New Mexico Example

Let’s go through the tree clump grouping process for a single stand flown on the Santa Fe National Forest in November 2024.

Setup

# library
library(tidyverse) # the tidyverse
library(viridis) # viridis colors
library(scales) # work with number and plot scales
library(latex2exp) # math formulas text
# library(paletteer) # working with color palettes
# library(palettetown) # color palettes
# palettetown::pokedex(cb = 4)
library(harrypotter) # color palettes
# spatial analysis
library(terra) # raster
library(sf) # simple features
# visualization
library(kableExtra)
library(mapview) # interactive html maps
library(patchwork) # ! better align plots in grid
library(ggnewscale) # add another scale to ggplot
library(ggtext) # customize plot text colors
library(ggpubr) # alternative plot fn's (table)
# forest patches
# install.packages("pak")
# pak::pkg_install("mhahsler/dbscan", upgrade = T)
library(dbscan)

now, define a function to convert metric to imperial units

calc_imperial_units_fn <- function(df) {
  df %>% 
  # convert to imperial units
    dplyr::mutate(
      dplyr::across(
        .cols = tidyselect::ends_with("_cm")
        , ~ .x * 0.394
        , .names = "{.col}_in"
      )
      , dplyr::across(
        .cols = tidyselect::ends_with("_m")
        , ~ .x * 3.28
        , .names = "{.col}_ft"
      )
      , dplyr::across(
        .cols = tidyselect::ends_with("_m2_per_ha")
        , ~ .x * 4.359
        , .names = "{.col}_ftac"
      )
      , dplyr::across(
        .cols = tidyselect::ends_with("_per_ha") & !tidyselect::ends_with("_m2_per_ha")
        , ~ .x * 0.405
        , .names = "{.col}_ac"
      )
      , dplyr::across(
        .cols = tidyselect::ends_with("_area_ha")
        , ~ .x * 2.471
        , .names = "{.col}_ac"
      )
      , dplyr::across(
        .cols = tidyselect::ends_with("_m2")
        , ~ .x * 10.764
        , .names = "{.col}_ft2"
      )
    ) %>%
    dplyr::rename_with(
      .fn = function(x){dplyr::case_when(
        stringr::str_ends(x,"_cm_in") ~ stringr::str_replace(x,"_cm_in","_in")
        , stringr::str_ends(x,"_m_ft") ~ stringr::str_replace(x,"_m_ft","_ft")
        , stringr::str_ends(x,"_m2_per_ha_ftac") ~ stringr::str_replace(x,"_m2_per_ha_ftac","_ft2_per_ac")
        , stringr::str_ends(x,"_per_ha_ac") ~ stringr::str_replace(x,"_per_ha_ac","_per_ac")
        , stringr::str_ends(x,"_area_ha_ac") ~ stringr::str_replace(x,"_area_ha_ac","_area_ac")
        , stringr::str_ends(x,"_m2_ft2") ~ stringr::str_replace(x,"_m2_ft2","_ft2")
        , TRUE ~ x
      )}
    )
}

Load Data

load data

# where is the processed data from cloud2trees?
input_dir <- "../data/"
# tree top points
treetops_sf_with_dbh <- sf::st_read(paste0(input_dir, "/final_detected_tree_tops.gpkg"))
## Reading layer `final_detected_tree_tops' from data source 
##   `C:\Data\usfs\rio_chama_cflrp_202411\data\final_detected_tree_tops.gpkg' 
##   using driver `GPKG'
## Simple feature collection with 80693 features and 17 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 341046.3 ymin: 4006169 xmax: 342265.5 ymax: 4007337
## Projected CRS: NAD83 / UTM zone 13N
# what?
treetops_sf_with_dbh %>% dplyr::glimpse()
## Rows: 80,693
## Columns: 18
## $ treeID                    <chr> "1_341396.7_4007336.9", "2_341418.7_4007336.…
## $ tree_height_m             <dbl> 14.132, 8.233, 10.531, 2.682, 5.747, 3.252, …
## $ crown_area_m2             <dbl> 1.16, 1.68, 0.56, 0.16, 0.36, 0.16, 0.16, 0.…
## $ fia_est_dbh_cm            <dbl> 32.854644, 15.114067, 20.838229, 4.806052, 9…
## $ fia_est_dbh_cm_lower      <dbl> 13.800244, 6.434515, 8.645597, 1.980668, 4.0…
## $ fia_est_dbh_cm_upper      <dbl> 58.496900, 26.748053, 37.532606, 8.588588, 1…
## $ dbh_cm                    <dbl> 32.854644, 15.114067, 20.838229, 4.806052, 9…
## $ is_training_data          <lgl> FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FA…
## $ dbh_m                     <dbl> 0.32854644, 0.15114067, 0.20838229, 0.048060…
## $ radius_m                  <dbl> 0.16427322, 0.07557033, 0.10419115, 0.024030…
## $ basal_area_m2             <dbl> 0.084778049, 0.017941244, 0.034104486, 0.001…
## $ basal_area_ft2            <dbl> 0.91255092, 0.19311955, 0.36710069, 0.019527…
## $ ptcld_extracted_dbh_cm    <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
## $ ptcld_predicted_dbh_cm    <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
## $ comp_trees_per_ha         <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
## $ comp_relative_tree_height <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
## $ comp_dist_to_nearest_m    <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
## $ geom                      <POINT [m]> POINT (341396.7 4007337), POINT (34141…
# where?
treetops_sf_with_dbh %>% 
  dplyr::slice_sample(prop = 0.01) %>% 
  mapview::mapview()

create a fake stand boundary

# let's just zoom in on the center ha and call this our "stand"
stand_sf <- treetops_sf_with_dbh %>% 
  sf::st_bbox() %>% 
  sf::st_as_sfc() %>%
  sf::st_centroid() %>% 
  sf::st_buffer(dist = sqrt(20000)/2, endCapStyle = "SQUARE")
# area
sf::st_area(stand_sf)
## 20000 [m^2]

check this stand on the map

zoom out if the map looks blank

mapview::mapviewOptions(
    homebutton = FALSE
    , basemaps = c("Esri.WorldImagery","OpenStreetMap")
  )
mapview::mapview(stand_sf, alpha.regions = 0, color = "red", lwd = 3)

zoom out if the map looks blank

let’s make a base plot of the stand

plt_stand <- ggplot() +
  geom_sf(data = stand_sf, color = "black", fill = NA) +
  scale_x_continuous(expand = c(0, 0)) +
  scale_y_continuous(expand = c(0, 0)) +
  labs(
    x = ""
    , y = ""
  ) +
  theme_void() +
  theme(
    legend.position = "top" # c(0.5,1)
    , legend.direction = "horizontal"
    , legend.margin = margin(0,0,0,0)
    , legend.text = element_text(size = 8)
    , legend.title = element_text(size = 8)
    , legend.key = element_rect(fill = NA, color = NA)
    # , plot.title = ggtext::element_markdown(size = 10, hjust = 0.5)
    , plot.title = element_text(size = 10, hjust = 0.5, face = "bold")
    , plot.subtitle = element_text(size = 8, hjust = 0.5, face = "italic")
  )

let’s define our DBH limit for classifying trees as “overstory” and we’ll define our distance in meters to separate tree clump groupings

tree_clump_dist_m <- 6
ostory_dbh_cm <- 6*2.54 # = inches*2.54

plot the stand with some trees

plt_stand +
  geom_sf(
    data = treetops_sf_with_dbh %>% 
      sf::st_intersection(stand_sf %>% sf::st_buffer(10)) %>% 
      # overstory
      dplyr::mutate(overstory = dbh_cm>=ostory_dbh_cm)
    , mapping = aes(color = overstory)
  )

ggsave("../data/01_tree_tops.jpg", width = 7, height = 7)

Define Clumps

create a general function to apply dbscan::dbscan to point data of class sf

With respect to clump size groupings, Churchill et al. (2016) note that:

Proportions for clump sizes should be lumped into four or five bins for operational simplicity. We use 4 or 5 bins (Fig 5): individual trees, small clumps (2-4 trees), medium clumps (5-9 trees), and large clumps (10-20+ trees). Note that when instructed to leave a large clump (e.g. 10-20 trees), marking crews often have difficulty leaving the upper end of the size range (e.g. an 18, 19, or 20 tree clump). Thus adding a fifth bin for “super clumps” may be necessary (e.g. 15-20 trees or 20-25+ trees), especially if the upper size range of clumps is desired. (p. 12-13)

In their prescription worksheet Excel file Churchill et al. (2016) have the following group sizes: 1, 2-4 (Mean clump size: 3), 5-9 (7), 10-15 (12), 16-25 (20). To be able to cut clumps down to achieve the desired largest “16-25” group, we’ll add a “super” group with >25 trees. This same principle can be applied if you desire a clump size group larger than this (i.e. 25-35), whereby a group above this should be added so that the resulting prescription cuts the “super” group to ensure the proper sizing of the desired largest group.

# function to clump data that are sf points
st_clump_points <- function(
  x # point data of class `sf` 
  , clump_dist_m = 6 # size (radius) of the epsilon neighborhood = maximum distance between points to add to clump
  , clump_breaks = c(0,1,4,9,15,25,Inf) # where to break the clump size groups
  , clump_labels = c("Individual","2-4 trees","5-9 trees","10-15 trees","16-25 trees",">25 trees") # what to name the clump size groups
) {
  # get points as x,y
    point_clusters = x %>% 
      dplyr::mutate(
        X = sf::st_coordinates(.)[,1] %>% as.numeric()
        , Y = sf::st_coordinates(.)[,2] %>% as.numeric()
      )
    #############################################################################
    ##### Identify clusters in each stem map plot                           #####
    #############################################################################
    ### Place trees into clusters using an inter-tree distance of 6 m
    my_dbscan_temp =  point_clusters %>% 
      sf::st_drop_geometry() %>% 
      dplyr::select(X,Y) %>% 
      dbscan::dbscan(eps = clump_dist_m, minPts = 2)
    
    # my_dbscan_temp %>% str()
    
    ### append cluster ID to tree points
    point_clusters$dbscan_cluster = my_dbscan_temp$cluster
    # point_clusters$cluster %>% summary()
    # point_clusters %>% sf::st_drop_geometry() %>% dplyr::count(cluster) %>% dplyr::arrange(desc(n)) %>% dplyr::slice_head(n=11)
  
  ### cluster metrics
  point_clusters = point_clusters %>% 
    dplyr::group_by(dbscan_cluster) %>% 
    dplyr::mutate(
      # unique dbscan_cluster for individuals
      clump_id = dplyr::case_when(
        dbscan_cluster == 0 ~ max(my_dbscan_temp$cluster)+dplyr::row_number()
        , T ~ dbscan_cluster
      ) %>% 
      factor()
    ) %>% 
    dplyr::group_by(clump_id) %>% 
    dplyr::mutate(
      dbscan_cluster = factor(dbscan_cluster)
      , clump_n_trees = dplyr::n()
      , clump_n_trees_grp = cut(
          clump_n_trees
          , breaks = clump_breaks
          , labels = clump_labels
        ) %>% 
        factor(
          ordered = T
          , levels = clump_labels
        )
    ) %>% 
    dplyr::ungroup() %>% 
    dplyr::mutate(tree_clump_dist_m = clump_dist_m)
  # return
  return(point_clusters)
}

function specific to this project to filter the tree points attached to harvest bounds and apply the st_clump_points function

# create function to pass a unit id and return list of trees with clump groupings
get_tree_clumps = function(
    tree_clump_dist_m=6
    , ostory_ht_m = as.numeric(NA)
    , ostory_dbh_cm = as.numeric(NA)
){
  # check ostory definition
  if(is.na(as.numeric(ostory_dbh_cm)) & is.na(as.numeric(ostory_ht_m))){
    warning("`ostory_dbh_cm` and `ostory_ht_m` are not set...using `ostory_dbh_cm` = 12.7")
    ostory_dbh_cm = 5*2.54
    # filter data
    ttops_temp = treetops_sf_with_dbh %>% sf::st_intersection(stand_sf) %>% 
      dplyr::filter(
        dbh_cm>=as.numeric(ostory_dbh_cm)
      )
  }else if(!is.na(as.numeric(ostory_dbh_cm))){
    # filter data
    ttops_temp = treetops_sf_with_dbh %>% sf::st_intersection(stand_sf) %>% 
      dplyr::filter(
        dbh_cm>=as.numeric(ostory_dbh_cm)
      )
  }else{
    # filter data
    ttops_temp = treetops_sf_with_dbh %>% sf::st_intersection(stand_sf) %>% 
      dplyr::filter(
        tree_height_m>=as.numeric(ostory_ht_m)
      )
  }
  
  ### Place trees into clusters using an inter-tree distance of 6 m
  ttops_temp = st_clump_points(
    x = ttops_temp
    , clump_dist_m = tree_clump_dist_m
  )
  
  # add distance to nearest within clump
  ttops_temp =
    ttops_temp %>% 
    dplyr::group_by(clump_id) %>%
    tidyr::nest() %>% 
    dplyr::mutate(
      distance_clump_nn_m = purrr::map(data, function(x){
        # get index of nearest neighbor
        i = sf::st_nearest_feature(x)
        # get dist
        d = sf::st_distance(x, x[i,], by_element=TRUE) %>% as.numeric()
        return(d)
      })
    ) %>% 
    tidyr::unnest(cols = c(data, distance_clump_nn_m)) %>% 
    sf::st_set_geometry("geom") %>% # set it cuz it got lost 
    dplyr::ungroup() %>% 
    dplyr::mutate(
      tree_clump_dist_m = tree_clump_dist_m
      , clump_id_duplicate = clump_id # can use this even after nesting data by clump_id
      # , ostory_ht_m = ifelse(is.na(ostory_ht_m), as.numeric(NA), as.numeric(ostory_ht_m))
      # , ostory_dbh_cm = ifelse(is.na(ostory_dbh_cm), as.numeric(NA), as.numeric(ostory_dbh_cm))
    )
  # return
  return(ttops_temp)
}

call that function

# call it
ttops_temp = get_tree_clumps(
  tree_clump_dist_m = 6
  , ostory_dbh_cm = ostory_dbh_cm # 6 in dbh
)
## Warning: attribute variables are assumed to be spatially constant throughout
## all geometries

what?

ttops_temp %>% dplyr::glimpse()
## Rows: 769
## Columns: 27
## $ clump_id                  <fct> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,…
## $ treeID                    <chr> "49736_341606.7_4006823.3", "49892_341603.5_…
## $ tree_height_m             <dbl> 20.767, 9.622, 15.808, 12.870, 11.495, 10.35…
## $ crown_area_m2             <dbl> 32.72, 1.52, 6.28, 4.28, 3.52, 1.68, 8.48, 2…
## $ fia_est_dbh_cm            <dbl> 61.52282, 18.49889, 38.74657, 28.38951, 23.8…
## $ fia_est_dbh_cm_lower      <dbl> 25.860036, 7.815217, 16.136208, 12.057785, 9…
## $ fia_est_dbh_cm_upper      <dbl> 109.37700, 32.93805, 68.86084, 51.30199, 42.…
## $ dbh_cm                    <dbl> 61.52282, 18.49889, 38.74657, 28.38951, 23.8…
## $ is_training_data          <lgl> FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FA…
## $ dbh_m                     <dbl> 0.6152282, 0.1849889, 0.3874657, 0.2838951, …
## $ radius_m                  <dbl> 0.30761410, 0.09249444, 0.19373285, 0.141947…
## $ basal_area_m2             <dbl> 0.29727772, 0.02687702, 0.11791157, 0.063300…
## $ basal_area_ft2            <dbl> 3.1998974, 0.2893042, 1.2692001, 0.6813641, …
## $ ptcld_extracted_dbh_cm    <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
## $ ptcld_predicted_dbh_cm    <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
## $ comp_trees_per_ha         <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
## $ comp_relative_tree_height <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
## $ comp_dist_to_nearest_m    <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
## $ geom                      <POINT [m]> POINT (341606.7 4006823), POINT (34160…
## $ X                         <dbl> 341606.7, 341603.5, 341609.9, 341616.9, 3416…
## $ Y                         <dbl> 4006823, 4006822, 4006822, 4006822, 4006821,…
## $ dbscan_cluster            <fct> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,…
## $ clump_n_trees             <int> 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, …
## $ clump_n_trees_grp         <ord> >25 trees, >25 trees, >25 trees, >25 trees, …
## $ tree_clump_dist_m         <dbl> 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6,…
## $ distance_clump_nn_m       <dbl> 3.352611, 2.154066, 3.352611, 2.630589, 2.72…
## $ clump_id_duplicate        <fct> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,…

Clump Polygons and Metrics

Churchill et al. (2016) provide instructions for implementing the clump identification (Plotkin et al. 2002) in ArcGIS:

Use the Buffer tool (in the Proximity toolset within the Analysis toolbox) to create a buffer of distance d/2, one half the inter-tree distance, around each point. This quantity d/2 is meant to approximate the crown radius of a “typical” overstory tree. Set the Dissolve Type option to ALL, which dissolves overlapping buffers, creating a reduced set of spatially non-overlapping polygons stored as a multipart polygon feature…Sanchez Meador et al. (2011) provide some useful examples of how clump attributes can be summarized…The method described here can be modified to use measured or modeled crown radii for each tree in place of d/2 (p.36)

# create function to pass a return from st_clump_points/get_tree_clumps and create clump polygons with summary stats
get_clump_summary = function(dta){
  # get tree_clump_dist_m
  tree_clump_dist_m = min(dta$tree_clump_dist_m, na.rm = T)
  # create clump polys and summary
  clump_polys_temp = 
    dta %>% 
      dplyr::ungroup() %>% 
      sf::st_set_geometry("geometry") %>% 
      sf::st_buffer(tree_clump_dist_m/2) %>% 
      dplyr::group_by(clump_id, dbscan_cluster, clump_n_trees_grp) %>%
      dplyr::summarise(
        # union buffered tree points
        geometry = sf::st_union(geometry)
        # summary metrics
        , n_trees = dplyr::n_distinct(treeID)
        , mean_dbh_cm = mean(dbh_cm, na.rm = T)
        , mean_tree_height_m = mean(tree_height_m, na.rm = T)
        , loreys_height_m = sum(basal_area_m2*tree_height_m, na.rm = T) / sum(basal_area_m2, na.rm = T)
        , basal_area_m2 = sum(basal_area_m2, na.rm = T)
        , sum_dbh_cm_sq = sum(dbh_cm^2, na.rm = T)
        , .groups = "drop_last"
      ) %>%
      dplyr::ungroup() %>% 
      sf::st_make_valid() %>% 
      dplyr::mutate(
        clump_area_ha = sf::st_area(geometry) %>% as.numeric() %>% `/`(10000)
        , trees_per_ha = (n_trees/clump_area_ha)
        , basal_area_m2_per_ha = (basal_area_m2/clump_area_ha)
        , pct_stand_basal_area = basal_area_m2/sum(basal_area_m2)
        , pct_stand_n_trees = n_trees/sum(n_trees)
        , qmd_cm = sqrt(sum_dbh_cm_sq/n_trees)
      ) %>%
      dplyr::select(-c(sum_dbh_cm_sq)) %>% 
    # convert to imperial units
      calc_imperial_units_fn() %>% 
      dplyr::mutate(tree_clump_dist_m = tree_clump_dist_m)
  # calculate distance between clumps
  clump_polys_temp = clump_polys_temp %>% 
    dplyr::mutate(
      nearest = sf::st_nearest_feature(clump_polys_temp)
      , distance_nearest_clump_m = sf::st_distance(
          clump_polys_temp
          , clump_polys_temp[nearest,]
          , by_element=TRUE
        ) %>% 
        as.numeric()
    ) %>% 
    dplyr::select(-c(nearest))
  # return
  return(clump_polys_temp)
}

get the clump summary data

# get it
# get_clump_summary(
#   dta = get_tree_clumps(tree_clump_dist_m = tree_clump_dist_m)
# )
clump_polys_temp = get_clump_summary(ttops_temp)
# what?
clump_polys_temp %>% dplyr::glimpse()
## Rows: 39
## Columns: 25
## $ clump_id                 <fct> 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 3…
## $ dbscan_cluster           <fct> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 2, 3, …
## $ clump_n_trees_grp        <ord> Individual, Individual, Individual, Individua…
## $ geometry                 <GEOMETRY [m]> POLYGON ((341593.9 4006821,..., POLY…
## $ n_trees                  <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 31, 83, 1…
## $ mean_dbh_cm              <dbl> 96.40704, 35.15035, 94.58870, 97.47311, 66.14…
## $ mean_tree_height_m       <dbl> 26.72900, 14.89200, 26.50900, 26.99600, 21.65…
## $ loreys_height_m          <dbl> 26.729000, 14.892000, 26.509001, 26.996000, 2…
## $ basal_area_m2            <dbl> 0.72997403, 0.09703962, 0.70269748, 0.7462073…
## $ clump_area_ha            <dbl> 0.002826142, 0.002826142, 0.002826142, 0.0028…
## $ trees_per_ha             <dbl> 353.8393, 353.8393, 353.8393, 353.8393, 353.8…
## $ basal_area_m2_per_ha     <dbl> 258.29350, 34.33643, 248.64199, 264.03750, 12…
## $ pct_stand_basal_area     <dbl> 0.0071688845, 0.0009530008, 0.0069010087, 0.0…
## $ pct_stand_n_trees        <dbl> 0.00130039, 0.00130039, 0.00130039, 0.0013003…
## $ qmd_cm                   <dbl> 96.40704, 35.15035, 94.58870, 97.47311, 66.14…
## $ mean_dbh_in              <dbl> 37.984375, 13.849236, 37.267948, 38.404406, 2…
## $ qmd_in                   <dbl> 37.984375, 13.849236, 37.267948, 38.404406, 2…
## $ mean_tree_height_ft      <dbl> 87.67112, 48.84576, 86.94952, 88.54688, 71.03…
## $ loreys_height_ft         <dbl> 87.67112, 48.84576, 86.94952, 88.54688, 71.03…
## $ basal_area_ft2_per_ac    <dbl> 1125.90139, 149.67250, 1083.83043, 1150.93946…
## $ trees_per_ac             <dbl> 143.3049, 143.3049, 143.3049, 143.3049, 143.3…
## $ clump_area_ac            <dbl> 0.006983396, 0.006983396, 0.006983396, 0.0069…
## $ basal_area_ft2           <dbl> 7.8574404, 1.0445344, 7.5638357, 8.0321761, 3…
## $ tree_clump_dist_m        <dbl> 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, …
## $ distance_nearest_clump_m <dbl> 2.12332457, 0.04846820, 0.01763506, 0.8018667…
# do these numbers match
identical(
  # clump 
  nrow(clump_polys_temp)
  # clumps in tree list data
  , ttops_temp %>% dplyr::distinct(clump_id) %>% nrow()
)
## [1] TRUE

plot the clump group size

plt_stand +
  geom_sf(
    data = clump_polys_temp
    , mapping = aes(fill = clump_n_trees_grp)
    , color = NA
  ) +
  harrypotter::scale_fill_hp_d("always", name = "clump size", drop = F) +
  labs(fill = "") +
  guides(size = "none", fill = guide_legend(override.aes = list(size = 5)))

ggsave("../data/02_tree_clumps.jpg", width = 7, height = 7)

Clump Spacing

See Churchill et al. (2016) Figure 4 (p.10) and Matonis and Binkley (2018) who “calculated coverage of mosaic-meadows (percentage of stand > 6 m from overstory trees)” (p. 124)

Since we already buffered the tree points to approximate the crown radius, we’ll continue to use our \(d/2\) where \(d\) is maximum distance between trees for determining tree clumps and is meant to approximate the crown radius of a “typical” overstory tree

# create function to pass a return from get_clump_summary() and get a distance raster
get_clump_dist_rast = function(dta){
  # get tree_clump_dist_m
  tree_clump_dist_m = min(dta$tree_clump_dist_m, na.rm = T)
  # rasterize the clump polygons and then calculate distance between clumps as raster
  dist_rast = 
    terra::rasterize(
        x = dta %>% terra::vect()
        , y = dta %>% 
          terra::vect() %>% 
          terra::rast(res = 0.2)
      ) %>% 
      terra::distance() %>% 
      # crop it to stand extent
      terra::crop(
        stand_sf %>% 
        terra::vect()
      ) %>% 
      terra::mask(
        stand_sf %>% 
        terra::vect()
      )
  ######### part 2
  # now create openings vector data
  openings_vect = 
    dist_rast %>% 
      terra::classify(rcl = c(tree_clump_dist_m/2,Inf), others = NA, include.lowest = T) %>% 
      terra::as.polygons(na.rm = T) %>% 
      sf::st_as_sf() %>% 
      sf::st_cast("POLYGON") %>% 
      dplyr::mutate(layer = dplyr::row_number()) %>% 
      dplyr::mutate(
        openining_area_m2 = sf::st_area(geometry) %>% as.numeric()
        , tree_clump_dist_m = tree_clump_dist_m
      )
  
  # return
  return(list(dist_rast = dist_rast, openings_vect = openings_vect))
}
# get it
dist_rast_temp = get_clump_dist_rast(clump_polys_temp)

what is it?

dist_rast_temp
## $dist_rast
## class       : SpatRaster 
## dimensions  : 708, 708, 1  (nrow, ncol, nlyr)
## resolution  : 0.2, 0.2  (x, y)
## extent      : 341585.1, 341726.7, 4006682, 4006824  (xmin, xmax, ymin, ymax)
## coord. ref. : NAD83 / UTM zone 13N (EPSG:26913) 
## source(s)   : memory
## name        :    layer 
## min value   :  0.00000 
## max value   : 13.64844 
## 
## $openings_vect
## Simple feature collection with 52 features and 3 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: 341585.1 ymin: 4006682 xmax: 341726.7 ymax: 4006824
## Projected CRS: NAD83 / UTM zone 13N
## First 10 features:
##     layer                       geometry openining_area_m2 tree_clump_dist_m
## 1       1 POLYGON ((341585.1 4006824,...              0.08                 6
## 1.1     2 POLYGON ((341596.5 4006824,...              2.88                 6
## 1.2     3 POLYGON ((341660.3 4006824,...              8.32                 6
## 1.3     4 POLYGON ((341622.5 4006824,...              7.08                 6
## 1.4     5 POLYGON ((341693.1 4006824,...             18.68                 6
## 1.5     6 POLYGON ((341676.1 4006818,...              2.48                 6
## 1.6     7 POLYGON ((341695.5 4006816,...              0.20                 6
## 1.7     8 POLYGON ((341662.5 4006802,...              0.04                 6
## 1.8     9 POLYGON ((341645.7 4006804,...              3.40                 6
## 1.9    10 POLYGON ((341662.1 4006801,...              0.04                 6

just view the distance raster

terra::plot(dist_rast_temp$dist_rast)

plot the distance raster and openings vector data we just got with overlaid tree clumps and tree points

plt_fnl_temp = 
  ggplot() + 
    # distance
    geom_tile(
      data = dist_rast_temp$dist_rast %>% terra::aggregate(2, cores = 4) %>% as.data.frame(xy = T) %>% rename(f=3)
      , mapping = aes(x=x, y=y, fill = f)
    ) +
    harrypotter::scale_fill_hp(
      "mischief"
      , na.value = "transparent"
      , name = "distance to\nnearest tree (m)"
      , limits = c(0,23)
    ) +
    # openings
    geom_sf(data = dist_rast_temp$openings_vect, mapping = aes(color = openining_area_m2), fill = NA) +
    scale_color_gradient(
      low = "gray77", high = "gray11"
      , labels = scales::comma_format(accuracy = 1)
      , name = latex2exp::TeX("opening\narea ($\\m^2$)")
      , limits = c(0,7500)
    ) +
    # clumps
    ggnewscale::new_scale_fill() +
    geom_sf(
      data = clump_polys_temp
      , mapping = aes(fill = clump_n_trees_grp)
      , color = NA
    ) +
    harrypotter::scale_fill_hp_d("always", name = "clump size", drop = F) +
    # tree points
    geom_sf(data = ttops_temp, color = "gray88", shape = ".") +
    theme_void() +
    theme(
      plot.subtitle = element_text(size = 10, hjust = 0.5, face = "bold")
      , legend.title=element_text(size=8), legend.text=element_text(size = 7)
    )
# plot
plt_fnl_temp

# save it
ggsave("../data/03_clumps_dist.jpg", width = 7, height = 7)

highlight the openings

plt_open_temp = 
  ggplot() + 
      # clumps
      geom_sf(
        data = clump_polys_temp
        , mapping = aes(fill = clump_n_trees_grp)
        , color = NA
      ) +
      harrypotter::scale_fill_hp_d("always", name = "clump size", drop = F) + 
      # openings
      ggnewscale::new_scale_fill() +
      geom_sf(data = dist_rast_temp$openings_vect, mapping = aes(fill = openining_area_m2), color = NA) +
      scale_fill_gradient(
        low = "gray77", high = "gray11"
        , labels = scales::comma_format(accuracy = 1)
        , name = latex2exp::TeX("opening\narea ($\\m^2$)")
        , limits = c(0,7500)
      ) +
      # tree points
      geom_sf(data = ttops_temp, color = "gray88", shape = ".") +
      theme_void() +
      theme(plot.subtitle = element_text(size = 10, hjust = 0.5, face = "bold"), legend.title=element_text(size=8), legend.text=element_text(size = 7))
# plot
plt_open_temp

ggsave("../data/04_clumps_openings.jpg", width = 7, height = 7)

combine them?

plt_fnl_temp + (plt_open_temp + theme(legend.position = "none")) &
  theme(plot.subtitle = element_text(size = 10, hjust = 0.5, face = "bold"), legend.title=element_text(size=8), legend.text=element_text(size = 7))

ggsave("../data/05_clumps_combined.jpg", width = 10, height = 7)

Clump Metrics

create a function to summarize by number of tree clump grouping variable

# create a function to summarize by number of tree clump grouping
get_clump_n_trees_grp_summary = function(trees, clumps){
  # get area of harvest unit
  #...will use this area in the area calculations such that...
  #...TPA = trees in a certain group size across the whole stand area
  harvest_area_ha = stand_sf %>% 
    sf::st_as_sf() %>% 
    dplyr::mutate(harvest_area_ha = sf::st_area(.) %>% as.numeric() %>% `/`(10000)) %>% 
    dplyr::pull(harvest_area_ha) %>% 
    .[1]
  # collapse and calculate silv metrics
  dta = trees %>% 
    sf::st_drop_geometry() %>% 
    dplyr::mutate(stand_area_ha = harvest_area_ha) %>% 
    dplyr::group_by(stand_area_ha,clump_n_trees_grp) %>%
    dplyr::summarise(
      # summary metrics
      n_trees = dplyr::n_distinct(treeID)
      , mean_dbh_cm = mean(dbh_cm, na.rm = T)
      , mean_tree_height_m = mean(tree_height_m, na.rm = T)
      , loreys_height_m = sum(basal_area_m2*tree_height_m, na.rm = T) / sum(basal_area_m2, na.rm = T)
      , basal_area_m2 = sum(basal_area_m2, na.rm = T)
      , sum_dbh_cm_sq = sum(dbh_cm^2, na.rm = T)
      , .groups = "drop_last"
    ) %>%
    dplyr::ungroup() %>% 
    # attach clump area
    dplyr::left_join(
      clumps %>% 
        sf::st_drop_geometry() %>% 
        dplyr::group_by(clump_n_trees_grp) %>%
        dplyr::summarise(
          clump_area_ha = sum(clump_area_ha)
          , stand_n_clumps = dplyr::n()
          , .groups = "drop_last"
        ) %>% 
        dplyr::ungroup()
      , by = dplyr::join_by(clump_n_trees_grp)
    ) %>% 
    dplyr::mutate(
      trees_per_ha = (n_trees/stand_area_ha) # (n_trees/clump_area_ha) ... this was not right
      , basal_area_m2_per_ha = (basal_area_m2/stand_area_ha) # (basal_area_m2/clump_area_ha) ... this was not right
      , qmd_cm = sqrt(sum_dbh_cm_sq/n_trees)
      # stand calcs
      , stand_trees_per_ha = sum(n_trees)/stand_area_ha
      , stand_basal_area_m2 = sum(basal_area_m2)
      , stand_basal_area_m2_per_ha = sum(basal_area_m2)/stand_area_ha
      , pct_stand_basal_area = basal_area_m2/stand_basal_area_m2
      , pct_stand_n_trees = n_trees/sum(n_trees)
      , stand_qmd_cm = sqrt(sum(trees$dbh_cm^2, na.rm = T)/sum(n_trees))
    ) %>%
    dplyr::select(-c(sum_dbh_cm_sq)) %>% 
    # convert to imperial units
    calc_imperial_units_fn()
  # return
  return(dta)
}

call it

# call it
clump_n_trees_grp_summary_temp = get_clump_n_trees_grp_summary(
  trees = get_tree_clumps(tree_clump_dist_m = tree_clump_dist_m, ostory_dbh_cm = ostory_dbh_cm)
  , clumps = get_clump_summary(
    get_tree_clumps(tree_clump_dist_m = tree_clump_dist_m, ostory_dbh_cm = ostory_dbh_cm)
  )
) 
# what?
clump_n_trees_grp_summary_temp %>% dplyr::glimpse()
## Rows: 6
## Columns: 31
## $ stand_area_ha               <dbl> 2, 2, 2, 2, 2, 2
## $ clump_n_trees_grp           <ord> Individual, 2-4 trees, 5-9 trees, 10-15 tr…
## $ n_trees                     <int> 12, 30, 31, 11, 40, 645
## $ mean_dbh_cm                 <dbl> 61.66965, 45.29162, 34.25379, 29.24668, 46…
## $ mean_tree_height_m          <dbl> 20.18100, 16.11127, 13.89961, 12.77318, 16…
## $ loreys_height_m             <dbl> 24.45865, 24.82122, 20.04323, 15.45925, 23…
## $ basal_area_m2               <dbl> 4.2887219, 7.2452159, 3.8558698, 0.8625523…
## $ clump_area_ha               <dbl> 0.03391370, 0.06613346, 0.05884907, 0.0159…
## $ stand_n_clumps              <int> 12, 13, 5, 1, 2, 6
## $ trees_per_ha                <dbl> 6.0, 15.0, 15.5, 5.5, 20.0, 322.5
## $ basal_area_m2_per_ha        <dbl> 2.1443610, 3.6226079, 1.9279349, 0.4312761…
## $ qmd_cm                      <dbl> 67.45721, 55.45237, 39.79563, 31.59741, 54…
## $ stand_trees_per_ha          <dbl> 384.5, 384.5, 384.5, 384.5, 384.5, 384.5
## $ stand_basal_area_m2         <dbl> 101.8253, 101.8253, 101.8253, 101.8253, 10…
## $ stand_basal_area_m2_per_ha  <dbl> 50.91266, 50.91266, 50.91266, 50.91266, 50…
## $ pct_stand_basal_area        <dbl> 0.042118420, 0.071153376, 0.037867491, 0.0…
## $ pct_stand_n_trees           <dbl> 0.01560468, 0.03901170, 0.04031209, 0.0143…
## $ stand_qmd_cm                <dbl> 41.06008, 41.06008, 41.06008, 41.06008, 41…
## $ mean_dbh_in                 <dbl> 24.29784, 17.84490, 13.49599, 11.52319, 18…
## $ qmd_in                      <dbl> 26.57814, 21.84823, 15.67948, 12.44938, 21…
## $ stand_qmd_in                <dbl> 16.17767, 16.17767, 16.17767, 16.17767, 16…
## $ mean_tree_height_ft         <dbl> 66.19368, 52.84495, 45.59073, 41.89604, 54…
## $ loreys_height_ft            <dbl> 80.22437, 81.41361, 65.74180, 50.70633, 76…
## $ basal_area_ft2_per_ac       <dbl> 9.347269, 15.790948, 8.403868, 1.879933, 2…
## $ stand_basal_area_ft2_per_ac <dbl> 221.9283, 221.9283, 221.9283, 221.9283, 22…
## $ trees_per_ac                <dbl> 2.4300, 6.0750, 6.2775, 2.2275, 8.1000, 13…
## $ stand_trees_per_ac          <dbl> 155.7225, 155.7225, 155.7225, 155.7225, 15…
## $ stand_area_ac               <dbl> 4.942, 4.942, 4.942, 4.942, 4.942, 4.942
## $ clump_area_ac               <dbl> 0.08380075, 0.16341579, 0.14541606, 0.0393…
## $ basal_area_ft2              <dbl> 46.163803, 77.987504, 41.504582, 9.284513,…
## $ stand_basal_area_ft2        <dbl> 1096.048, 1096.048, 1096.048, 1096.048, 10…

summary table

# table it
clump_n_trees_grp_summary_temp %>% 
  dplyr::select(
    clump_n_trees_grp, n_trees
    , mean_dbh_in
    , qmd_in
    , mean_tree_height_ft
    , loreys_height_ft
    , trees_per_ac
    , basal_area_ft2_per_ac, pct_stand_basal_area, pct_stand_n_trees
  ) %>% 
  dplyr::mutate(
    dplyr::across(
      .cols = c(pct_stand_basal_area, pct_stand_n_trees) 
      , .fns = ~ scales::percent(.x, accuracy = 1)
    )
  ) %>% 
  kableExtra::kbl(
    digits = 1
    , escape = F
    , caption = paste0("Overstory tree clump summary")
    , col.names = c(
      "", "trees"
      , "mean<br>DBH (in)"
      , "QMD (in)"
      , "mean<br>Ht. (ft)"
      , "Loreys<br>Ht. (ft)"
      , "TPA"
      , "BA<br>ft<sup>2</sup> ac<sup>-1</sup>"
      , "%<br>stand BA"
      , "%<br>stand trees"
      )
  ) %>% 
  kableExtra::kable_styling()
Overstory tree clump summary
trees mean
DBH (in)
QMD (in) mean
Ht. (ft)
Loreys
Ht. (ft)
TPA BA
ft2 ac-1
%
stand BA
%
stand trees
Individual 12 24.3 26.6 66.2 80.2 2.4 9.3 4% 2%
2-4 trees 30 17.8 21.8 52.8 81.4 6.1 15.8 7% 4%
5-9 trees 31 13.5 15.7 45.6 65.7 6.3 8.4 4% 4%
10-15 trees 11 11.5 12.4 41.9 50.7 2.2 1.9 1% 1%
16-25 trees 40 18.3 21.4 54.6 76.7 8.1 20.2 9% 5%
>25 trees 645 13.6 15.3 46.2 61.6 130.6 166.3 75% 84%

ICO Implementation

Churchill et al. (2016) describe the full process for implementing the ICO approach in The ICO Approach to Quantifying and Restoring Forest Spatial Pattern: Implementation Guide in which the authors lay out the prescription development process:

  1. Identify skips and other special treatment areas
  2. Consider the need for openings
  3. Determine the stand average density target
  4. Determine the appropriate distance to define clumps
  5. Obtain targets for clump proportions
  6. Select target clump proportions for your stand
  7. Generate clump targets for the whole unit
  8. Combine clump and opening targets with leave tree criteria into marking guidelines

The objective here is to: 1) provide the manager with the current conditions (completed above); 2) take the “targets” as set by the manager (steps 3, 5, 6, 7); 3) create the prescription with the leave tree marking.

Let’s implement this prescription development process with our UAS tree list

3. Determine the stand average density target

Step 3 in Churchill et al. (2016):

An average BA, TPA, or SDI target for the stand should be selected that is appropriate for the species, structure, site conditions, and management objectives. Expected mortality from prescribed fire should be factored in. Stand average targets can come from historical reference stands, plant association based stocking guides, density management tools, or a combination of both (see Franklin et al. (2013) for a full discussion of setting density targets). In dry forests, the number and size of old trees must be accounted in setting the density target. To use the ICO method, the target must be converted to TPA (see Table 1). A lower diameter cutoff also needs to be specified for the TPA target. This should be the lower limit in the contract or cutting guidelines given to the marking crew or contractor. (p.11)

this is what Table 1 looks like with TPA values are derived from the formula:

\[ TPA = \frac{BA}{QMD^{2} \times 0.005454} \]

# function to get tpa from ba and qmd
get_tpa = function(ba_ft2_ac, qmd_in){
  tpa = round(ba_ft2_ac/((qmd_in^2)*0.005454))
  return(tpa)
} 
# table it
tidyr::crossing(
    ba = seq(40,200,20)
    , qmd = seq(8,20,2)
  ) %>% 
  dplyr::mutate(
    tpa = get_tpa(ba,qmd)
  ) %>% 
  tidyr::pivot_wider(names_from = ba, values_from = tpa) %>% 
  dplyr::mutate(l = "QMD (in)") %>% 
  dplyr::relocate(l) %>% 
  kableExtra::kbl(
    col.names = c(".","", seq(40,200,20))
    , escape = F
    , caption = "Basal Area and QMD to TPA conversion chart"
  ) %>% 
  kableExtra::add_header_above(
    c("","", "Basal Area (ft2/ac)"=length(seq(40,200,20)))
  ) %>%
  kableExtra::kable_styling() %>% 
  kableExtra::column_spec(1:2, bold = T) %>%
  kableExtra::collapse_rows(columns = 1, valign = "middle")
Basal Area and QMD to TPA conversion chart
Basal Area (ft2/ac)
. 40 60 80 100 120 140 160 180 200
QMD (in) 8 115 172 229 286 344 401 458 516 573
10 73 110 147 183 220 257 293 330 367
12 51 76 102 127 153 178 204 229 255
14 37 56 75 94 112 131 150 168 187
16 29 43 57 72 86 100 115 129 143
18 23 34 45 57 68 79 91 102 113
20 18 28 37 46 55 64 73 83 92

Current Stand Conditions

For determining targets, the silviculturist needs to know the current conditions. Provide the current stand conditions based on the UAS tree list for the selected stand that are required to set the targets:

  • Current BA
  • Current QMD
  • Current proportion of trees by clump size
clump_n_trees_grp_summary_temp %>% 
  dplyr::select(clump_n_trees_grp, pct_stand_n_trees) %>% 
  dplyr::mutate(
    pct_stand_n_trees = scales::percent(pct_stand_n_trees,accuracy = 1)
  ) %>% 
  tidyr::pivot_wider(names_from = clump_n_trees_grp, values_from = pct_stand_n_trees) %>% 
  kableExtra::kable(
    caption = paste0(
      "Current stand BA (ft2/ac): "
      , clump_n_trees_grp_summary_temp$stand_basal_area_ft2_per_ac[1] %>% scales::number(accuracy = 0.1)
      , "<br>Current stand QMD (in): "
      , clump_n_trees_grp_summary_temp$stand_qmd_in[1] %>% scales::number(accuracy = 0.1)
      , "<br>Current stand TPA: "
      , clump_n_trees_grp_summary_temp$stand_trees_per_ac[1] %>% scales::number(accuracy = 1)
    )
    , escape = F
    , digits = 1
  ) %>% 
  kableExtra::kable_styling() %>% 
  kableExtra::footnote(general = "values are the percent of trees in each clump size")
Current stand BA (ft2/ac): 221.9
Current stand QMD (in): 16.2
Current stand TPA: 156
Individual 2-4 trees 5-9 trees 10-15 trees 16-25 trees >25 trees
2% 4% 4% 1% 5% 84%
Note:
values are the percent of trees in each clump size

5. Obtain targets for clump proportions

Step 5 in Churchill et al. (2016):

ICO prescriptions are based on a target proportion of trees in different sized clumps within a stand. Proportions are just the percentage of trees, or TPA, that are in different sized clumps. Basal area proportions can be used, but we have found TPA targets to be more straightforward to use. Ideally, a table summarizing clump proportions for a range of reference conditions in your area is available (Table 2). If not, instructions for developing one are provided in section VI. (p.12)

Section VI of Churchill et al. (2016) notes that

reference spatial information may already be available and summarized in a way that it can be directly incorporated into ICO prescriptions. Such data exist and have been published for areas in Arizona (Abella and Denton 2009, Sánchez Meador et al. 2011), the eastern Washington Cascades (Churchill et al. 2013), the northern Rockies (Larson et al. 2012), and the Sierra Nevada (Lydersen et al. 2013). Reference datasets for using ICO in other forest types, such as coastal Douglas-fir or Pacific silver fir, also exist (Larson and Churchill 2008). (p.28)

Table 2 is:

6. Select target clump proportions for your stand

Now set the desired BA, QMD, and proportion of trees in each clump size:

########################################################################################
########################################################################################
# desired BA, QMD, and proportion of trees in each clump size
########################################################################################
########################################################################################
# desired BA
  target_ba = 85 # cannot be > current BA
# desired QMD
  target_qmd = 18
# desired proportion (%) of trees in each clump size
  # !cannot be create larger proportion of ">25 trees" clump as this would require adding trees...
  # c("Individual", "2-4 trees",    "5-9 trees",    "10-15 trees","16-25 trees",">25 trees")
  # c(.18, .33, .24, .10, .15)
  target_pcts = c(30, 25, 20, 10, 15, 0)
########################################################################################
########################################################################################
# desired BA, QMD, and proportion of trees in each clump size
########################################################################################
########################################################################################

Check set up and define data with targets based on Churchill et al. (2016) and their prescription worksheet Excel file to help develop ICO prescriptions

get_target_check_prescription = function(
  clump_n_trees_grp_summary_dta
  , target_ba = as.numeric(NA)
  , target_qmd = as.numeric(NA)
  , target_pcts = as.numeric(NA)
){
  if(
    is.na(target_ba) | is.na(target_qmd) | max(is.na(target_pcts))==1
  ){
    stop("must set all of the function parameters:\n`target_ba`, `target_qmd`, and `target_pcts`")
  }
  #############################################
  # check target BA and TPA
  #############################################
  if(as.numeric(target_ba)>clump_n_trees_grp_summary_dta$stand_basal_area_ft2_per_ac[1]){
    stop(
      "target BA in `target_ba` of "
      , round(as.numeric(target_ba),1), " is greater than current BA of "
      , clump_n_trees_grp_summary_dta$stand_basal_area_ft2_per_ac[1] %>% round(1)
    )
  }
  if(
    get_tpa(target_ba, target_qmd)>clump_n_trees_grp_summary_dta$stand_trees_per_ac[1]
  ){
    stop(
      "target TPA in of "
      , round(as.numeric(get_tpa(target_ba, target_qmd)),1), " is greater than current TPA of "
      , clump_n_trees_grp_summary_dta$stand_trees_per_ac[1] %>% round(1)
      , "\n adjust `target_ba` and/or `target_qmd` to get valid TPA"
    )
  }
  #############################################
  # define data with current and target
  # ... this is "smart" in that percentages are adj based on:
  # ... 0) are there missing targets?
  # ... ... if < 6 numbers provided in `target_pcts` then the largest tree groups get targets of 0
  # ... 1) do targets sum to 1? 
  # ... ... if not trees are distributed proportionally based on targets provided and trees available
  # ... 2) is target in largest clump size > current conditions?
  # ... ... if yes, target is set to current condition
  # ... 3) is target listed in clump size > current largest clump with trees?
  # ... ... if yes, target for largest clump size is shifted to current largest clump with trees
  #############################################
    target_data = 
      # create data for joining if missing clump groups
      dplyr::tibble(
        stand_area_ac = rep(clump_n_trees_grp_summary_dta$stand_area_ac[1],6)
        , clump_n_trees_grp = factor(
          c(1:6)
          , labels = c("Individual", "2-4 trees",   "5-9 trees",    "10-15 trees","16-25 trees",">25 trees")
          , ordered = T
        )
        , mean_clump_n_trees = c(1,3,7,12,20,30)
        , min_clump_n_trees = c(1,2,5,10,16,26)
        , max_clump_n_trees = c(1,4,9,15,25,99999)
      ) %>% 
      dplyr::left_join(
        clump_n_trees_grp_summary_dta %>% 
          dplyr::ungroup() %>% 
          dplyr::select(clump_n_trees_grp, pct_stand_n_trees, stand_n_clumps)
      ) %>% 
      dplyr::mutate(
        pct_stand_n_trees = dplyr::coalesce(pct_stand_n_trees,0)
        , stand_n_clumps = dplyr::coalesce(stand_n_clumps,0)
      ) %>% 
      # add targets
      dplyr::bind_cols(
        pct_stand_n_trees_target = c(as.numeric(target_pcts), rep(0,6))[1:6] # pad target with 0's
      ) %>% 
      # adjust target based on difference from 1
      dplyr::mutate(
        pct_stand_n_trees_target = pct_stand_n_trees_target*(1/sum(pct_stand_n_trees_target))
        # largest clump size with trees
        , largest_w_trees = max(ifelse(dplyr::coalesce(pct_stand_n_trees)>0,clump_n_trees_grp,NA),na.rm = T)
        , largest_w_trees_target = max(ifelse(dplyr::coalesce(pct_stand_n_trees_target)>0,clump_n_trees_grp,NA),na.rm = T)
      ) %>% 
      # move target for largest clump size to the largest current clump size 
      dplyr::mutate(
        pct_stand_n_trees_target = dplyr::case_when(
          as.numeric(clump_n_trees_grp)==largest_w_trees &
            largest_w_trees_target>largest_w_trees ~ max(
              ifelse(as.numeric(clump_n_trees_grp)==largest_w_trees_target,pct_stand_n_trees_target,0)
            )
          , T ~ pct_stand_n_trees_target
        )
      ) %>% 
      # adjust target based on current conditions
      dplyr::mutate(
        pct_stand_n_trees_target = dplyr::case_when(
          as.numeric(clump_n_trees_grp)>largest_w_trees &
            pct_stand_n_trees_target > 0 ~ 0
          , as.numeric(clump_n_trees_grp)==largest_w_trees &
            pct_stand_n_trees_target > pct_stand_n_trees ~ pct_stand_n_trees
          , T ~ pct_stand_n_trees_target
        )
      ) %>% 
      # finally, re-scale again based on adjustments
      dplyr::mutate(
        pct_stand_n_trees_target = dplyr::case_when(
          as.numeric(clump_n_trees_grp)==largest_w_trees ~ pct_stand_n_trees_target
          , T ~ pct_stand_n_trees_target * (
            # pct remaining to scale to
            (1-max(ifelse(as.numeric(clump_n_trees_grp)==largest_w_trees,pct_stand_n_trees_target,0))) /
            # current pct remaining total allocated
            sum(
              ifelse(as.numeric(clump_n_trees_grp)!=largest_w_trees,pct_stand_n_trees_target,0))
            )
        )
      ) %>% 
      # add other targets
      dplyr::rename(pct_stand_n_trees_current = pct_stand_n_trees) %>% 
      dplyr::mutate(
        stand_trees_per_ac_current = clump_n_trees_grp_summary_dta$stand_trees_per_ac[1]
        , stand_trees_per_ac_target = dplyr::coalesce(get_tpa(target_ba, target_qmd),0)
        , trees_per_acre_current = stand_trees_per_ac_current*pct_stand_n_trees_current
        , trees_per_acre_target = stand_trees_per_ac_target*pct_stand_n_trees_target
        , clumps_per_acre_current = trees_per_acre_current/mean_clump_n_trees
        , clumps_per_acre_target = trees_per_acre_target/mean_clump_n_trees
        , stand_n_clumps_current = stand_n_clumps
        , stand_n_clumps_target = (clumps_per_acre_target*stand_area_ac) %>% round(0)
      ) %>% 
      dplyr::select(-c(tidyselect::starts_with("largest_w_trees"), stand_n_clumps))
  # ????  
  # target_data %>% glimpse()
  
    # issue warning about targets
    if(min(target_data$pct_stand_n_trees_target == c(as.numeric(target_pcts), rep(0,6))[1:6])==0){
      warning(
        "proportion of trees in each clump size target `target_pcts` adjusted!!!"
        , "\nfrom : ", paste(round(target_pcts,2),collapse = ",")
        , "\nto : ", paste(round(target_data$pct_stand_n_trees_target,2),collapse = ",")
      )
    }
  # return
  return(target_data)
}
# call it
target_data_temp = get_target_check_prescription(
  clump_n_trees_grp_summary_temp
  , target_ba = target_ba
  , target_qmd = target_qmd
  , target_pcts = target_pcts
)
## Warning in get_target_check_prescription(clump_n_trees_grp_summary_temp, : proportion of trees in each clump size target `target_pcts` adjusted!!!
## from : 30,25,20,10,15,0
## to : 0.3,0.25,0.2,0.1,0.15,0
# what?
target_data_temp %>% dplyr::glimpse()
## Rows: 6
## Columns: 15
## $ stand_area_ac              <dbl> 4.942, 4.942, 4.942, 4.942, 4.942, 4.942
## $ clump_n_trees_grp          <ord> Individual, 2-4 trees, 5-9 trees, 10-15 tre…
## $ mean_clump_n_trees         <dbl> 1, 3, 7, 12, 20, 30
## $ min_clump_n_trees          <dbl> 1, 2, 5, 10, 16, 26
## $ max_clump_n_trees          <dbl> 1, 4, 9, 15, 25, 99999
## $ pct_stand_n_trees_current  <dbl> 0.01560468, 0.03901170, 0.04031209, 0.01430…
## $ pct_stand_n_trees_target   <dbl> 0.30, 0.25, 0.20, 0.10, 0.15, 0.00
## $ stand_trees_per_ac_current <dbl> 155.7225, 155.7225, 155.7225, 155.7225, 155…
## $ stand_trees_per_ac_target  <dbl> 48, 48, 48, 48, 48, 48
## $ trees_per_acre_current     <dbl> 2.4300, 6.0750, 6.2775, 2.2275, 8.1000, 130…
## $ trees_per_acre_target      <dbl> 14.4, 12.0, 9.6, 4.8, 7.2, 0.0
## $ clumps_per_acre_current    <dbl> 2.4300000, 2.0250000, 0.8967857, 0.1856250,…
## $ clumps_per_acre_target     <dbl> 14.400000, 4.000000, 1.371429, 0.400000, 0.…
## $ stand_n_clumps_current     <dbl> 12, 13, 5, 1, 2, 6
## $ stand_n_clumps_target      <dbl> 71, 20, 7, 2, 2, 0

current vs target

target_data_temp %>% 
  dplyr::select(
    clump_n_trees_grp
    , tidyselect::starts_with("pct_stand_n_trees")
    , tidyselect::starts_with("trees_per_acre_")
  ) %>% 
  dplyr::mutate(
    dplyr::across(
      tidyselect::starts_with("pct_stand_n_trees")
      , ~ scales::percent(.x,accuracy = 1)
    )
  ) %>% 
  kableExtra::kable(
    caption = paste0(
      "Current stand BA (ft2/ac): "
      , clump_n_trees_grp_summary_temp$stand_basal_area_ft2_per_ac[1] %>% scales::number(accuracy = 0.1)
      , "<br>Current stand QMD (in): "
      , clump_n_trees_grp_summary_temp$stand_qmd_in[1] %>% scales::number(accuracy = 0.1)
      , "<br>Current stand TPA: "
      , clump_n_trees_grp_summary_temp$stand_trees_per_ac[1] %>% scales::number(accuracy = 1)
    )
    , escape = F
    , digits = 1
    , col.names = c(
      "", rep(c("current","target"),2)
    )
  ) %>% 
  kableExtra::kable_styling() %>% 
  kableExtra::add_header_above(
    c(" "=1,"% Trees"=2, "TPA"=2)
  )
Current stand BA (ft2/ac): 221.9
Current stand QMD (in): 16.2
Current stand TPA: 156
% Trees
TPA
current target current target
Individual 2% 30% 2.4 14.4
2-4 trees 4% 25% 6.1 12.0
5-9 trees 4% 20% 6.3 9.6
10-15 trees 1% 10% 2.2 4.8
16-25 trees 5% 15% 8.1 7.2
>25 trees 84% 0% 130.6 0.0

8. Combine clump and opening targets with leave tree criteria into marking guidelines

Use our UAS tree list to generate the prescription:

  • start with the largest clump size currently with trees
  • cut trees to the next largest clump size until desired # clumps is reached
  • repeat with each successive clump size through to individual tree selection
  • if possible, cut in same clump until desired proportions are reached to minimize machine time

Function to cut clumps

put the entire process outlined immediately above into a function to use with the tree list to cut each clump and then we’ll select the desired proportion of clumps

define intermediate functions

define intermediate functions

##############################################
# working with sf LINESTRINGS
##############################################
  # first two functions borrowed from https://github.com/metafor-ulaval/ALSroads/blob/main/R/line_tools.R
  ########
  # Get heading of both ends of a line
  ########
    st_ends_heading <- function(line){
      M <- sf::st_coordinates(line)
      i <- c(2, nrow(M) - 1)
      j <- c(1, -1)
      
      headings <- mapply(i, j, FUN = function(i, j) {
        Ax = M[i-j,1]
        Ay = M[i-j,2]
        Bx = M[i,1]
        By = M[i,2]
        atan2(Ay-By, Ax-Bx)*180/pi
      })
      names(headings) <- c("head", "tail")
      return(headings)
    }
  ########
  # extend the line on both ends
  ########
    st_extend_line <- function(line, distance, end = "BOTH"){
      if (!(end %in% c("BOTH", "HEAD", "TAIL")) | length(end) != 1) stop("'end' must be 'BOTH', 'HEAD' or 'TAIL'")
    
      M <- sf::st_coordinates(line)[,-3]
      keep <- !(end == c("TAIL", "HEAD"))
      
      ends <- c(1, nrow(M))[keep]
      headings <- st_ends_heading(line)[keep] / 180 * pi
      distances <- if (length(distance) == 1) rep(distance, 2) else distance[1:2]
      
      M[ends, 1:2] <- M[ends, 1:2] + distances[keep] * c(cos(headings), sin(headings))
      newline <- sf::st_linestring(M)
    
      # If input is sfc_LINESTRING and not sfg_LINESTRING
      if (is.list(line)) newline <- sf::st_sfc(newline, crs = sf::st_crs(line))
      
      return(newline)
    }
  ########
  # pass an sf dataframe of points and return a line between the farthest points
  ########
    st_points_to_line <- function(pts, line_ext=0) {
      if(max(class(ttops_temp) %in% c("sf"))!=1){
        stop("must provide an object of class `sf`")
      }
      # find farthest distance between points
      dist_temp = sf::st_distance(pts)
      
      # get the points
      f_pts_temp = 
        pts %>% 
          dplyr::ungroup() %>% 
          dplyr::slice(
            # get the farthest points from distance matrix
            which(dist_temp == max(dist_temp), arr.ind = TRUE)[1,]
          )
      
      # draw a line between the farthest two points
      f_line_temp = f_pts_temp %>% 
          # convert to linestring
          dplyr::ungroup() %>% 
          dplyr::summarise(n=dplyr::n()) %>%
          sf::st_cast("LINESTRING") %>%
          dplyr::ungroup() %>% 
          dplyr::select(-c(n))
      
      # and apply the line extension
        farthest_line = st_extend_line(f_line_temp, distance = line_ext)
      # return
        return(farthest_line)
    }
    # st_points_to_line(ttops_temp %>% dplyr::slice_head(prop = .1), line_ext = 6) %>% 
    #   ggplot() + geom_sf() + geom_sf(data = ttops_temp %>% dplyr::slice_head(prop = .1)) + theme_void()

  ########
  # Function to calculate Euclidean distance between 2 points
  ########
    st_euclidean_distance <- function(p1,p2) {
        return(sqrt((p2[1] - p1[1])^2 + (p2[2] - p1[2])^2))
    }

  ########
  # return a line perpendicular to current line
  ########
  ### https://stackoverflow.com/questions/56771058/perpendicular-lines-at-regular-intervals-along-lines-with-multiple-nodes
    # Function to calculate 2 points on a line perpendicular to another defined by 2 points p1,p2
    # For point at interval, which can be a proportion of the segment length, or a constant
    st_perp_line <- function(interval=0.5, my_line, proportion=TRUE) {
      # get end points of line
      p1 = my_line %>% sf::st_cast("POINT") %>% sf::st_coordinates() %>% .[1,]
      p2 = my_line %>% sf::st_cast("POINT") %>% sf::st_coordinates() %>% .[2,]
      # get length of line to return equal length line
      line_len = sf::st_length(my_line) %>% as.numeric() %>% `/`(2)      
      # get crs of line
      my_crs = sf::st_crs(my_line)
      
      # Calculate x and y distances
      x_len <- p2[1] - p1[1]
      y_len <- p2[2] - p1[2]
      
      # If proportion calculate reference point from tot_length
      if (proportion) {
          point <- c(p1[1]+x_len*interval,p1[2]+y_len*interval)
      }
      # Else use the constant value
      else {
          tot_len <- st_euclidean_distance(p1,p2)
          point <- c(p1[1]+x_len/tot_len*interval,p1[2]+y_len/tot_len*interval)
      }    
      
      # Calculate the x and y distances from reference point to point on line line_len distance away    
      ref_len <- st_euclidean_distance(point,p2)
      xn_len <- (line_len / ref_len) * (p2[1] - point[1])
      yn_len <- (line_len / ref_len) * (p2[2] - point[2])
      
      # fix for identical 
      if(identical(point,p2) & x_len>y_len){ # this works for horizontal line
        xn_len <- line_len/2
        yn_len <- 0
      }else if(identical(point,p2) & x_len<y_len){ # this works for vertical line
        xn_len <- 0
        yn_len <- line_len/2
      }
      
      # Invert the x and y lengths and add/subtract from the refrence point
      # ref_points <- rbind(point,c(point[1] + yn_len,point[2] - xn_len),c(point[1] - yn_len,point[2] + xn_len))
      ref_points <- rbind(c(point[1] + yn_len,point[2] - xn_len),c(point[1] - yn_len,point[2] + xn_len))
      
      # use the reference points to return a line
      return(
        ref_points %>% 
          dplyr::as_tibble() %>% 
          dplyr::rename_with(tolower) %>% 
          sf::st_as_sf(coords = c("x", "y"), crs = my_crs, remove = F) %>%
          dplyr::summarise(n=dplyr::n()) %>%
          sf::st_cast("LINESTRING") %>%
          dplyr::ungroup() %>% 
          dplyr::select(-c(n))
      )
    }
function to pass a clump

generate data using clumping functions above, pass that data, return tree list with cut/keep flag based on the target clump group size

# function to cut a single clump
cut_clump_fn <- function(
      c # clump id from need_cut_trees data 
      , need_cut_trees # data with clumps already defined returned by st_clump_points()
      , tgt = tgt # select one level of input passed to st_clump_points() ...
          # ... clump_labels = c("Individual","2-4 trees","5-9 trees","10-15 trees","16-25 trees",">25 trees")
){
  # use only the current clump
  curr_dta = need_cut_trees %>% 
    dplyr::filter(clump_id == c)
  
  if(nrow(curr_dta)<1){
    stop("cannot find data with the clump_id == `c` in `need_cut_trees` data")
  }
  
  # get distance used to create clumps in st_clump_points()
  dist_temp = curr_dta$tree_clump_dist_m[1]
  
  # do we even need to cut?
  reclump = st_clump_points(x = curr_dta, clump_dist_m = dist_temp)
  if(
    (
      reclump %>% 
        dplyr::pull(clump_id) %>% 
        unique() %>% 
        length()
    ) != 1
  ){
    stop("this is not a clump...send the `need_cut_trees` data through the st_clump_points function again")
  }
  # do we even need to cut?
  if(
    unique(reclump$clump_n_trees_grp) == tgt
  ){
    return(
      reclump %>% 
        sf::st_drop_geometry() %>% 
        dplyr::select(treeID, is_keep_tree)
    )
  }
  
  # create clump polygon
  clumps = get_clump_summary(curr_dta)
  
  # get the farthest line between points
  f_line_temp = st_points_to_line(curr_dta, line_ext = dist_temp)
  
  # get perpendicular lines in dataset which we can iterate over to make cuts
  perp_line_sf_temp = 
    # for every 1 m along line length, get a new perp line
    seq(
        from = 0
        , to = sf::st_length(f_line_temp) %>% as.numeric() %>% floor()
        , by = 1
      ) %>% 
    purrr::map(
      st_perp_line
      , my_line = f_line_temp
      , proportion = F
    ) %>% 
    dplyr::bind_rows() %>% 
    dplyr::mutate(line_n = dplyr::row_number())
  
  # find intersection of lines with the polygon and add length of intersection to perp line data
  perp_line_sf_temp = perp_line_sf_temp %>% 
    dplyr::inner_join(
      # intersect and calc len
      perp_line_sf_temp %>% 
        sf::st_intersection(
          clumps %>% 
            dplyr::ungroup() %>% 
            dplyr::select(clump_id) %>% 
            dplyr::filter(clump_id == c)
        ) %>% 
        dplyr::mutate(len_m = sf::st_length(geometry) %>% as.numeric()) %>% 
        sf::st_drop_geometry()
      , by = "line_n"
    )
  
  # make cuts at the points where there is the least overlap with the clump polygon
    # list of potential line cut + tree combinations
    # aggregate the mean length of the intersecting cut lines to the tree level
    # prioritize trees for removal that have smallest length
    cut_tree_lines_temp = 
      curr_dta %>% 
        sf::st_buffer(dist_temp/2) %>% 
        sf::st_join(perp_line_sf_temp) %>% 
        sf::st_drop_geometry() %>% 
        dplyr::group_by(treeID) %>% 
        dplyr::summarise(len_m = mean(len_m, na.rm = T)) %>% 
        dplyr::ungroup() %>% 
        dplyr::arrange(len_m, treeID) %>% 
        dplyr::mutate(n = dplyr::row_number())
  ###############################################################
  # while
  ###############################################################
  # cut until the desired clump size is achieved based on these cut lines
  while_temp = 1
  i = 1
  while(while_temp==1) {
    # cut trees
    cut_trees_temp = cut_tree_lines_temp %>% 
      dplyr::slice(1:i) %>% 
      dplyr::distinct(treeID)
    
    # get the remaining trees not cut
    trees_remain_temp = curr_dta %>% 
      dplyr::anti_join(cut_trees_temp, by = "treeID")
    
    # ensure that there are trees
    if(nrow(trees_remain_temp)==0){
      if(best_desired_grps_n_temp==0){
        # there are no more possible cuts :/
        # ... so we're going to add trees until the desired clump size is reached
        # start with biggest tree and keep adding trees until the desired clump size is reached
          start_tree = curr_dta %>% 
            dplyr::arrange(desc(dbh_cm)) %>% 
            dplyr::filter(dplyr::row_number() == 1)
          
          increment_m = 0.5
          k = 1
          keep_trees = dplyr::tibble(treeID = character(0))
          while(nrow(keep_trees)==0){
            # create a polygon to intersect with tree points and keep trees until number of trees met
            poly_keep = start_tree %>% 
              sf::st_buffer( (dist_temp/2) + increment_m*k ) %>% 
              dplyr::mutate(dummy = 1) %>% 
              dplyr::select(dummy)
            # intersect and clump
            i_trees = curr_dta %>% 
              sf::st_intersection(poly_keep) %>% 
              st_clump_points(clump_dist_m = dist_temp)
            i_trees %>% sf::st_drop_geometry() %>% dplyr::count(clump_n_trees_grp, clump_id)
            # keep only the desired clump
            keep_trees = i_trees %>% 
              sf::st_drop_geometry() %>% 
              dplyr::filter(
                clump_id == 
                  i_trees %>% 
                    dplyr::filter(clump_n_trees_grp == tgt) %>% 
                    dplyr::pull(clump_id) %>% 
                    .[1] %>% 
                    dplyr::coalesce("nope")
              ) %>% 
              dplyr::select(treeID)
            # increment
            k = k+1
          }
          # the remaining trees are cuts
          best_cuts = curr_dta %>% 
            sf::st_drop_geometry() %>% 
            dplyr::anti_join(keep_trees, by = "treeID") %>% 
            dplyr::select(treeID)
      }else if(
        best_desired_grps_n_temp>0 
      ){ 
        ###############################
        # add trees back in until gets worse
        ###############################
        j = 1
        while_add = 1
        while(while_add==1){
          # cut trees
          cut_trees_temp = best_cuts %>% 
            # add trees back in (i.e. remove them from the cut trees)
            dplyr::anti_join(
              cut_tree_lines_temp %>% 
                dplyr::slice(1:j) %>% 
                dplyr::distinct(treeID)
              , by = "treeID"
            )
          # get the remaining trees not cut
          trees_remain_temp = curr_dta %>% 
            dplyr::anti_join(cut_trees_temp, by = "treeID")
          
          # count the groups remaining after cuts
          desired_grps_n_temp = trees_remain_temp %>% 
            st_clump_points(clump_dist_m = dist_temp) %>% 
            sf::st_drop_geometry() %>% 
            # do we have group sizes we want?
            dplyr::filter(clump_n_trees_grp == tgt) %>% 
            nrow()
          
          if(desired_grps_n_temp>=best_desired_grps_n_temp){ # is this better than the best
            best_cuts = cut_trees_temp
            best_desired_grps_n_temp = desired_grps_n_temp
          }else{
            # stop it
            while_add = 0  
          }
          # increment
          j = j + 1
        } # while(while_add==1)
      } # best_desired_grps_n_temp>0 
      # done so stop the whole stop it
        while_temp = 0
    }else{ # if(nrow(trees_remain_temp)==0)
      # count the groups remaining after cuts
      desired_grps_n_temp = trees_remain_temp %>% 
        st_clump_points(clump_dist_m = dist_temp) %>% 
        sf::st_drop_geometry() %>% 
        # do we have group sizes we want?
        dplyr::filter(clump_n_trees_grp == tgt) %>% 
        nrow()
      ### store best cut list
      if(i==1){
        best_cuts = cut_trees_temp
        best_desired_grps_n_temp = desired_grps_n_temp
      }else if(desired_grps_n_temp>best_desired_grps_n_temp){ # is this better than the best
        best_cuts = cut_trees_temp
        best_desired_grps_n_temp = desired_grps_n_temp
      }else if(
        desired_grps_n_temp==best_desired_grps_n_temp
        & i!=nrow(cut_tree_lines_temp)
      ){
        best_cuts = best_cuts
        best_desired_grps_n_temp = best_desired_grps_n_temp
      }else if(
        best_desired_grps_n_temp>0 
        & desired_grps_n_temp<best_desired_grps_n_temp
      ){ # is this worse than the best which was successful?
        ###############################
        # add trees back in until gets worse
        ###############################
        j = 1
        while_add = 1
        while(while_add==1){
          # cut trees
          cut_trees_temp = best_cuts %>% 
            # add trees back in (i.e. remove them from the cut trees)
            dplyr::anti_join(
              cut_tree_lines_temp %>% 
                dplyr::slice(1:j) %>% 
                dplyr::distinct(treeID)
              , by = "treeID"
            )
          # get the remaining trees not cut
          trees_remain_temp = curr_dta %>% 
            dplyr::anti_join(cut_trees_temp, by = "treeID")
          
          # count the groups remaining after cuts
          desired_grps_n_temp = trees_remain_temp %>% 
            st_clump_points(clump_dist_m = dist_temp) %>% 
            sf::st_drop_geometry() %>% 
            # do we have group sizes we want?
            dplyr::filter(clump_n_trees_grp == tgt) %>% 
            nrow()
          
          if(desired_grps_n_temp>=best_desired_grps_n_temp){ # is this better than the best
            best_cuts = cut_trees_temp
            best_desired_grps_n_temp = desired_grps_n_temp
          }else{
            # stop it
            while_add = 0  
          }
          # increment
          j = j + 1
        } # while(while_add==1)
        # done so stop the whole stop it
        while_temp = 0
      }else if( i==nrow(cut_tree_lines_temp) ){ # is this the end?
        # stop it
        while_temp = 0
      }
      ### increment
      i = i+1
    } # else if(nrow(trees_remain_temp)==0)
  } # while(while_temp==1)
  
  # return it
  # return treelist with cut/keep
  # join to original data and pull
  d_temp = curr_dta %>%
    sf::st_drop_geometry() %>%
    dplyr::mutate(is_keep_tree = as.numeric(NA)) %>% 
    dplyr::select(-c(is_keep_tree)) %>% 
    dplyr::left_join(
      best_cuts %>% dplyr::mutate(is_keep_tree = 0)
      , by = dplyr::join_by("treeID")
    ) %>%
    dplyr::mutate(is_keep_tree = dplyr::coalesce(is_keep_tree, 1)) %>%
    dplyr::select(treeID, is_keep_tree)
  # returns treeID and keep tree flag data frame
  return(d_temp)
}
# # example
# # pass it a clump id in the data
# ttops_temp %>%
#   dplyr::filter(clump_n_trees_grp == "2-4 trees") %>%
#   dplyr::arrange(desc(clump_n_trees)) %>%
#   dplyr::slice(1) %>%
#   dplyr::pull(clump_id) %>%
#   cut_clump_fn(
#     need_cut_trees = ttops_temp
#     , tgt = "Individual"
#   )

there are some cases where the cut_clumps_fn may not return the desired tree clumps based on a single pass using the initial cut lines…iterate over the function until all clumps are in the desired clump size or smaller.

get_cut_clump <- function(
  c # clump id from need_cut_trees data 
  , need_cut_trees # data with clumps already defined returned by st_clump_points()
  , tgt = tgt # select one level of input passed to st_clump_points() ...
      # ... clump_labels = c("Individual","2-4 trees","5-9 trees","10-15 trees","16-25 trees",">25 trees")
) {
  # what if we need to keep cutting because could not find a solution based on initial cut lines?
  cuts_first = cut_clump_fn(
    c = c
    , need_cut_trees = need_cut_trees
    , tgt = tgt
  )
  # did we get the target?
  new_clumps = 
    need_cut_trees %>% 
      dplyr::inner_join(
        cuts_first %>% dplyr::filter(is_keep_tree == 1) %>% dplyr::select(treeID)
        , by = "treeID"
      ) %>%
      st_clump_points()
  # new_clumps %>% sf::st_drop_geometry() %>% dplyr::count(clump_n_trees_grp)
  
  while(
    (
      new_clumps %>% dplyr::filter(clump_n_trees_grp > tgt) %>% nrow()
    ) > 0
  ){
    # redo the cut clump
    # redo the cut clump
    cuts_again = new_clumps %>% 
      dplyr::filter(clump_n_trees_grp > tgt) %>%
      dplyr::pull(clump_id) %>% 
      unique() %>% 
      purrr::map(cut_clump_fn, need_cut_trees = new_clumps, tgt = tgt) %>% 
      dplyr::bind_rows() %>% 
      dplyr::rename(updt = is_keep_tree)
    
    # update the original cut data
    cuts_first = cuts_first %>%
      dplyr::left_join(
        cuts_again
        , by = "treeID"
      ) %>% 
      dplyr::mutate(is_keep_tree = dplyr::coalesce(updt, is_keep_tree)) %>% 
      dplyr::select(-c(updt))
    
    # reset the new clumps
    new_clumps = 
      need_cut_trees %>% 
        dplyr::inner_join(
          cuts_first %>% dplyr::filter(is_keep_tree == 1) %>% dplyr::select(treeID)
          , by = "treeID"
        ) %>%
        st_clump_points()
    # new_clumps %>% sf::st_drop_geometry() %>% dplyr::count(clump_n_trees_grp)
  } # end while
  return(cuts_first) 
}
function to pass a tree list

function to pass a whole tree list of sf point data

##############################################
# get cut keep tree flag
##############################################
get_keep_tree_flag <- function(
  x # x = sf point data
  , tgt # tgt = clump_n_trees_grp level defined in st_clump_points: ...
    ## ... clump_labels = c("Individual","2-4 trees","5-9 trees","10-15 trees","16-25 trees",">25 trees")
  , clump_dist_m = 6 # size (radius) of the epsilon neighborhood = maximum distance between points to add to clump
) {
  # MAKE CUTS INDEPENDENT OF CURRENT CLUMP GROUPING...
  # what if we pass a tree list with trees already cut? and thus, potentially different clump group sizes than is 
  # currently defined in clump_n_trees_grp?
  # 1) apply the clump grouping 
  # 2) only go through the cut alg for clumps that need cutting to the target
  # 3) append the trees in target clump or lower to the list at the end
  
  #############################
  # 1) apply the clump grouping 
  #############################
  new_clump_trees = x %>% 
    st_clump_points(clump_dist_m = clump_dist_m)
    
  # need cutting still
  need_cut_trees = new_clump_trees %>% 
    dplyr::filter(clump_n_trees_grp > tgt)
  
  #############################
  # 2) only go through the cut alg for clumps that need cutting to the target
  #############################
  # for each clump that still needs cutting, iterate over and return keep/cut flag
  compl_need_cut_trees = need_cut_trees %>% 
    dplyr::pull(clump_id) %>% 
    unique() %>% 
    purrr::map(get_cut_clump, need_cut_trees = need_cut_trees, tgt = tgt) %>%  # purrr:map fn
    dplyr::bind_rows()
  
  #############################
  # 4) append the trees in target clump or lower to the list at the end
  #############################
  compl_need_cut_trees = compl_need_cut_trees %>% 
    dplyr::bind_rows(
      new_clump_trees %>% 
        sf::st_drop_geometry() %>% 
        dplyr::filter(clump_n_trees_grp <= tgt) %>% 
        dplyr::mutate(is_keep_tree = 1) %>%
        dplyr::select(treeID, is_keep_tree)
    )
  
  # return
  return(
    # original data so that order is preserved
    x %>% 
      # add on trees that got cut with trees already good
      dplyr::inner_join(
        compl_need_cut_trees
        , by = "treeID"
      ) %>% 
      dplyr::pull(is_keep_tree)
  )
  
} # get_keep_tree_flag

Process to cut groups incrementally

cut groups incrementally to achieve target tree clump group size proportions

# must first have target data
# target_data_temp = get_target_check_prescription(
#   clump_n_trees_grp_summary_temp
#   , target_ba = target_ba
#   , target_qmd = target_qmd
#   , target_pcts = target_pcts
# )

# what is the smallest group with a target?
sm_grp_temp = target_data_temp %>% 
    dplyr::filter(pct_stand_n_trees_target>0) %>% 
    dplyr::pull(clump_n_trees_grp) %>% 
    min()

# start with the largest group in target data currently with trees
grp_temp = target_data_temp %>% 
  dplyr::arrange(desc(clump_n_trees_grp)) %>% 
  dplyr::filter(pct_stand_n_trees_current>0) %>% 
  dplyr::slice(1) %>% 
  dplyr::pull(clump_n_trees_grp)

message(
  paste("started cutting", grp_temp, "at", Sys.time())
)
# identify clumps to leave untouched
# clump_polys_temp = get_clump_summary(ttops_temp)
keep_clumps_temp = 
  clump_polys_temp %>% 
    sf::st_drop_geometry() %>% 
    #keep only trees in current group
    dplyr::inner_join(
      target_data_temp %>% 
        dplyr::filter(clump_n_trees_grp == grp_temp) %>% 
        dplyr::select(
          clump_n_trees_grp, mean_clump_n_trees, min_clump_n_trees, max_clump_n_trees
          , stand_n_clumps_target
        )
      , by = "clump_n_trees_grp"
    ) %>% 
    # keep only clumps that meet criteria
    dplyr::filter(
      n_trees >= min_clump_n_trees
      & n_trees <= max_clump_n_trees
    ) %>% 
    # keep clumps closest the mean number in target
    dplyr::mutate(
      pct_to_target = abs( (n_trees-mean_clump_n_trees)/mean_clump_n_trees )
    ) %>% 
    dplyr::arrange(pct_to_target, desc(qmd_cm), n_trees, clump_id) %>% 
    dplyr::filter(dplyr::row_number()<=stand_n_clumps_target) %>% 
    dplyr::select(clump_id)

# !!!!!!!!!!FIX: WHAT IF WE HAVE CLUMPS OF THIS SIZE BUT THE N_TREES>MAX_TREES AND NEED TO GET CLUMPS OF THIS SIZE?
#   ... UPDATE CUTTING ALG TO CUT CLUMP DOWN TO MIN-MAX TREE RANGE FOR CLUMPS WITH INF UPPER LIMIT

# start building tree list with keep/cut flag
# ttops_temp = get_tree_clumps(
#   my_suid = my_suid
#   , tree_clump_dist_m = tree_clump_dist_m
#   , ostory_dbh_cm = ostory_dbh_cm
# )
# build tree list 
keep_trees =
  ttops_temp %>% 
    dplyr::ungroup() %>% 
    sf::st_drop_geometry() %>% 
    dplyr::inner_join(keep_clumps_temp, by = "clump_id") %>% 
    dplyr::select(treeID) %>% 
    dplyr::mutate(is_keep_tree = 1)

# did we keep all of the clumps?
if( 
  nrow(keep_clumps_temp) == 
  ( clump_polys_temp %>% dplyr::filter(clump_n_trees_grp == grp_temp) %>% nrow() )
){
  remaining_trees = dplyr::tibble(
    treeID = character(0)
    , is_keep_tree = numeric(0)
  )
}else{
  # determine keep/cut for the remaining trees in that group type 
  remaining_trees = 
    ttops_temp %>%
      dplyr::ungroup() %>% 
      # keep only trees in current grp
      dplyr::filter(clump_n_trees_grp == grp_temp) %>% 
      # remove keep trees
      dplyr::anti_join(keep_trees, by = "treeID")
  
  # apply the cutting algorithm
  # get the flag
  remaining_trees$is_keep_tree = get_keep_tree_flag(
      x = remaining_trees
      , tgt = target_data_temp %>%  
        dplyr::ungroup() %>%
        dplyr::arrange(clump_n_trees_grp) %>%
        dplyr::mutate(l = dplyr::lag(clump_n_trees_grp)) %>%
        dplyr::filter(clump_n_trees_grp == grp_temp) %>%
        dplyr::pull(l)
      , clump_dist_m = tree_clump_dist_m
    ) 
  
  # select relevant columns
  remaining_trees = remaining_trees %>% 
    sf::st_drop_geometry() %>% 
    dplyr::ungroup() %>% 
    dplyr::select(treeID, is_keep_tree)
}
message(
  paste("done cutting", grp_temp, "at", Sys.time())
)
###############################################
# now process to go on to the next groups
###############################################
# get the next group
  grp_temp = target_data_temp %>%
    dplyr::filter(clump_n_trees_grp<grp_temp) %>% 
    dplyr::arrange(desc(clump_n_trees_grp)) %>% 
    dplyr::slice(1) %>% 
    dplyr::pull(clump_n_trees_grp)

while(grp_temp>=sm_grp_temp & grp_temp != "Individual"){
  message(
    paste("started cutting", grp_temp, "at", Sys.time())
  )
  # identify clumps to leave untouched
  # clump_polys_temp = get_clump_summary(ttops_temp)
  keep_clumps_temp =
    clump_polys_temp %>% 
      sf::st_drop_geometry() %>% 
      #keep only trees in current group
      dplyr::inner_join(
        target_data_temp %>% 
          dplyr::filter(clump_n_trees_grp == grp_temp) %>% 
          dplyr::select(
            clump_n_trees_grp, mean_clump_n_trees, min_clump_n_trees, max_clump_n_trees
            , stand_n_clumps_target
          )
        , by = "clump_n_trees_grp"
      ) %>% 
      # keep only clumps that meet criteria
      dplyr::filter(
        n_trees >= min_clump_n_trees
        & n_trees <= max_clump_n_trees
      ) %>% 
      # keep clumps closest the mean number in target group
      dplyr::mutate(
        pct_to_target = abs( (n_trees-mean_clump_n_trees)/mean_clump_n_trees )
      ) %>% 
      dplyr::arrange(pct_to_target, desc(qmd_cm), n_trees, clump_id) %>% 
      dplyr::filter(dplyr::row_number()<=stand_n_clumps_target) %>% 
      dplyr::select(clump_id)
  
  # add to tree list with keep/cut flag
  keep_trees = keep_trees %>% 
    dplyr::bind_rows(
      ttops_temp %>% 
        dplyr::ungroup() %>% 
        sf::st_drop_geometry() %>% 
        dplyr::inner_join(keep_clumps_temp, by = "clump_id") %>% 
        dplyr::select(treeID) %>% 
        dplyr::mutate(is_keep_tree = 1)
    )

  # check if the desired clump number was met and add the remaining trees from previous group if needed
  more_clumps_target = 0
  if(
    nrow(keep_clumps_temp) <
    ( 
      target_data_temp %>% 
        dplyr::filter(clump_n_trees_grp == grp_temp) %>% 
        dplyr::pull(stand_n_clumps_target) 
    )
  ){
    # how many more clumps are needed?
    more_clumps_target = 
      ( 
        target_data_temp %>% 
          dplyr::filter(clump_n_trees_grp == grp_temp) %>% 
          dplyr::pull(stand_n_clumps_target) 
      ) - nrow(keep_clumps_temp)
    
    # determine group size of remaining trees in the group prior that were not in a group selected and were not cut
    potential_trees = 
      ttops_temp %>% 
        dplyr::ungroup() %>% 
        dplyr::inner_join(
          remaining_trees %>% dplyr::filter(is_keep_tree == 1) %>% dplyr::select(treeID)
          , by = "treeID"
        ) %>% 
        st_clump_points(clump_dist_m = tree_clump_dist_m) %>% 
        # ggplot() + geom_sf(aes(fill = clump_n_trees_grp)) + theme_void()
        # get the original clump id to prioritize new clump groups in the same area
        dplyr::inner_join(
          ttops_temp %>% 
            sf::st_drop_geometry() %>% 
            dplyr::ungroup() %>% 
            dplyr::select(treeID, clump_id) %>% 
            dplyr::rename(orig_clump_id = clump_id)
          , by = "treeID"
        ) %>% 
        dplyr::group_by(orig_clump_id) %>% 
        dplyr::mutate(
          pct_desired_grp = 
            sum(ifelse(clump_n_trees_grp == grp_temp, 1, 0)) / dplyr::n()
        ) %>% 
        dplyr::ungroup() %>% 
        # keep only the current group
        dplyr::filter(clump_n_trees_grp == grp_temp)
    
    # pick trees from potential trees based on clump summary
    keep_trees = keep_trees %>% 
      dplyr::bind_rows(
        potential_trees %>% 
          sf::st_drop_geometry() %>% 
          # filter trees based on clumps needed
          dplyr::inner_join(
            get_clump_summary(potential_trees) %>% 
              sf::st_drop_geometry() %>% 
              # join with original clump id metrics to prioritize 
              # keeping clumps in the same area and minimize cutting time
              dplyr::left_join(
                potential_trees %>% 
                  sf::st_drop_geometry() %>% 
                  dplyr::group_by(clump_id, orig_clump_id) %>% 
                  dplyr::summarise(pct_desired_grp = max(pct_desired_grp))
                , by = "clump_id"
              ) %>% 
              # keep the number of clumps needed
              dplyr::arrange(desc(pct_desired_grp), orig_clump_id, desc(qmd_cm), n_trees, clump_id) %>% 
              dplyr::filter(dplyr::row_number()<=more_clumps_target) %>%
              dplyr::select(clump_id)
            , by = "clump_id"
          ) %>% 
          dplyr::select(treeID) %>% 
          dplyr::mutate(is_keep_tree = 1)
      )
  } # end if don't have enough clumps
  ####################################################
  # determine keep/cut for the remaining trees in that group type and higher groups 
  ####################################################
  remaining_trees =
    ttops_temp %>%
      dplyr::ungroup() %>% 
      # REMOVE TREES FROM PREVIOUS TREES REMAINING that got cut to make this group size
      dplyr::anti_join(
        remaining_trees %>% dplyr::filter(is_keep_tree == 0)
        , by = "treeID"
      ) %>%
      # keep only trees in current grp or prior grp
      dplyr::filter(clump_n_trees_grp >= grp_temp) %>% 
      # remove keep trees
      dplyr::anti_join(keep_trees, by = "treeID")
  
  if(nrow(remaining_trees) > 0){
    # apply the cutting algorithm
    # get the flag
    remaining_trees$is_keep_tree = get_keep_tree_flag(
        x = remaining_trees
        , tgt = target_data_temp %>%  
          dplyr::ungroup() %>%
          dplyr::arrange(clump_n_trees_grp) %>%
          dplyr::mutate(l = dplyr::lag(clump_n_trees_grp)) %>%
          dplyr::filter(clump_n_trees_grp == grp_temp) %>%
          dplyr::pull(l)
        , clump_dist_m = tree_clump_dist_m
      ) 
    
    # select relevant columns
    remaining_trees = remaining_trees %>% 
      sf::st_drop_geometry() %>% 
      dplyr::ungroup() %>% 
      dplyr::select(treeID, is_keep_tree)
  }
  
  # increment
    message(
      paste("done cutting", grp_temp, "at", Sys.time())
    )
  # get the next group
  grp_temp = target_data_temp %>%
    dplyr::filter(clump_n_trees_grp<grp_temp) %>% 
    dplyr::arrange(desc(clump_n_trees_grp)) %>% 
    dplyr::slice(1) %>% 
    dplyr::pull(clump_n_trees_grp) %>% 
    dplyr::coalesce("Individual")
}
# now individuals
if(sm_grp_temp == "Individual"){
  message(
    paste("started cutting", sm_grp_temp, "at", Sys.time())
  )
  potential_trees = 
    # original data
    ttops_temp %>% 
    dplyr::filter(clump_n_trees_grp == "Individual") %>% 
    # remaining trees
    dplyr::bind_rows(
      ttops_temp %>% 
        dplyr::ungroup() %>% 
        dplyr::inner_join(
          remaining_trees %>% dplyr::filter(is_keep_tree == 1) %>% dplyr::select(treeID)
          , by = "treeID"
        )
    ) %>% 
    # make sure that these are all individuals
    dplyr::mutate(is_orig = ifelse(clump_n_trees_grp=="Individual", 1, 0)) %>% 
    st_clump_points(clump_dist_m = tree_clump_dist_m) %>% 
    dplyr::filter(clump_n_trees_grp == "Individual")
  
    # pick trees from potential trees based on target
    keep_trees = keep_trees %>% 
      dplyr::select(treeID) %>% 
      dplyr::bind_rows(
        potential_trees %>% 
          sf::st_drop_geometry() %>% 
          # keep the number of clumps needed
          dplyr::arrange(desc(is_orig), desc(dbh_cm), desc(tree_height_m)) %>% 
          dplyr::filter(
            dplyr::row_number() <=
              target_data_temp %>% 
              dplyr::filter(clump_n_trees_grp == "Individual") %>% 
              dplyr::pull(stand_n_clumps_target)
          ) %>% 
          dplyr::select(treeID)
      ) %>% 
      dplyr::mutate(is_keep_tree = 1)
    
  message(
    paste("done cutting", sm_grp_temp, "at", Sys.time())
  )
}
# get the final prescription
prescription_trees = ttops_temp %>% 
  dplyr::ungroup() %>% 
  dplyr::left_join(
    keep_trees
    , by = "treeID"
  ) %>% 
  # tracking vars
  dplyr::mutate(
    # fill in keep tree flag
    is_keep_tree = dplyr::coalesce(is_keep_tree, 0)
    , orig_clump_n_trees_grp = clump_n_trees_grp
  ) %>% 
  dplyr::select(-c(clump_n_trees_grp))

# attach the new clumping to the trees
prescription_trees = prescription_trees %>% 
  dplyr::left_join(
    # reclump
    prescription_trees %>% 
      dplyr::filter(is_keep_tree == 1) %>% 
      st_clump_points(clump_dist_m = tree_clump_dist_m) %>% 
      sf::st_drop_geometry() %>% 
      dplyr::select(treeID, clump_n_trees_grp) %>% 
      dplyr::rename(new_clump_n_trees_grp = clump_n_trees_grp)
    , by = "treeID"
  ) %>% 
  dplyr::mutate(
    new_clump_n_trees_grp = forcats::fct_na_value_to_level(new_clump_n_trees_grp, level = "Cut tree")
  )

check the distribution of current vs prescription tree clump sizes

prescription_trees %>% 
  sf::st_drop_geometry() %>% 
  dplyr::count(
    orig_clump_n_trees_grp, new_clump_n_trees_grp
  ) %>% 
  tidyr::pivot_wider(names_from = orig_clump_n_trees_grp, values_from = n, values_fill = 0) %>% 
  dplyr::arrange(new_clump_n_trees_grp) %>% 
  dplyr::rename(`New` = new_clump_n_trees_grp) %>% 
  kableExtra::kbl() %>% 
  kableExtra::kable_styling() %>% 
  kableExtra::add_header_above(
    c(" "=1,"Old"=6)
  )
Old
New Individual 2-4 trees 5-9 trees 10-15 trees 16-25 trees >25 trees
Individual 12 0 0 0 0 52
2-4 trees 0 30 0 0 0 9
5-9 trees 0 0 31 0 0 25
10-15 trees 0 0 0 11 0 10
16-25 trees 0 0 0 0 40 17
Cut tree 0 0 0 0 0 532

get the distance raster and openings vector data

ttops_temp2 = prescription_trees %>% 
  dplyr::filter(is_keep_tree == 1) %>% 
  dplyr::select(treeID, dbh_cm, tree_height_m, basal_area_m2) %>% 
  st_clump_points(clump_dist_m = tree_clump_dist_m)
clump_polys_temp2 = get_clump_summary(ttops_temp2)
dist_rast_temp2 = get_clump_dist_rast(clump_polys_temp2)
# what are the openings
dist_rast_temp2$openings_vect$openining_area_m2 %>% summary()
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    0.04    0.08    0.98  380.88   22.97 7460.28
dist_rast_temp2$dist_rast %>% terra::aggregate(2, cores = 4) %>% terra::summary()
##      layer        
##  Min.   : 0.0000  
##  1st Qu.: 0.1707  
##  Median : 2.1610  
##  Mean   : 3.2247  
##  3rd Qu.: 4.9830  
##  Max.   :23.4648

plot the prescription

plt_fnl_temp2 = 
  ggplot() + 
    # distance
    geom_tile(
      data = dist_rast_temp2$dist_rast %>% terra::aggregate(2, cores = 4) %>% as.data.frame(xy = T) %>% rename(f=3)
      , mapping = aes(x=x, y=y, fill = f)
    ) +
    harrypotter::scale_fill_hp(
      "mischief"
      , na.value = "transparent"
      , name = "distance to\nnearest tree (m)"
      , limits = c(0,23)
    ) +
    # openings
    geom_sf(data = dist_rast_temp2$openings_vect, mapping = aes(color = openining_area_m2), fill = NA) +
    scale_color_gradient(
      low = "gray77", high = "gray11"
      , labels = scales::comma_format(accuracy = 1)
      , name = latex2exp::TeX("opening\narea ($\\m^2$)")
      , limits = c(0,7500)
    ) +
    # clumps
    ggnewscale::new_scale_fill() +
    geom_sf(
      data = clump_polys_temp2
      , mapping = aes(fill = clump_n_trees_grp)
      , color = NA
    ) +
    harrypotter::scale_fill_hp_d("always", name = "clump size", drop = F) + 
    # tree points
    geom_sf(data = ttops_temp2, color = "gray88", shape = ".") +
    theme_void() +
    theme(plot.subtitle = element_text(size = 10, hjust = 0.5, face = "bold"), legend.title=element_text(size=8), legend.text=element_text(size = 7))
# plot
plt_fnl_temp2

# save it
ggsave("../data/06_clumps_dist_presciption.jpg", width = 7, height = 7)

highlight the openings

plt_open_temp2 = 
  ggplot() + 
      # clumps
      geom_sf(
        data = clump_polys_temp2
        , mapping = aes(fill = clump_n_trees_grp)
        , color = NA
      ) +
      harrypotter::scale_fill_hp_d("always", name = "clump size", drop = F) + 
      # openings
      ggnewscale::new_scale_fill() +
      geom_sf(data = dist_rast_temp2$openings_vect, mapping = aes(fill = openining_area_m2), color = NA) +
      scale_fill_gradient(
        low = "gray77", high = "gray11"
        , labels = scales::comma_format(accuracy = 1)
        , name = latex2exp::TeX("opening\narea ($\\m^2$)")
        , limits = c(0,7500)
      ) +
      # tree points
      geom_sf(data = ttops_temp2, color = "gray88", shape = ".") +
      theme_void() +
      theme(plot.subtitle = element_text(size = 10, hjust = 0.5, face = "bold"), legend.title=element_text(size=8), legend.text=element_text(size = 7))
plt_open_temp2

# save it
ggsave("../data/07_clumps_openings_presciption.jpg", width = 7, height = 7)

combine plots with prescription

plt_fnl_temp2 + (plt_open_temp2 + theme(legend.position = "none")) &
  theme(plot.subtitle = element_text(size = 10, hjust = 0.5, face = "bold"), legend.title=element_text(size=8), legend.text=element_text(size = 7))

ggsave("../data/08_clumps_combined_prescription.jpg", width = 10, height = 7)

combine old vs new

(plt_fnl_temp + labs(subtitle = "Pre-Treatment")) +
  (plt_fnl_temp2  + labs(subtitle = "Post-Treatment") + theme(legend.position = "none")) &
  theme(plot.subtitle = element_text(size = 10, hjust = 0.5, face = "bold"), legend.title=element_text(size=8), legend.text=element_text(size = 7))

# save it
ggsave("../data/09_clumps_dist_compared.jpg", width = 10, height = 7)

combine old vs new openings highlight

(plt_open_temp + labs(subtitle = "Pre-Treatment")) +
  (plt_open_temp2  + labs(subtitle = "Post-Treatment") + theme(legend.position = "none")) &
  theme(plot.subtitle = element_text(size = 10, hjust = 0.5, face = "bold"), legend.title=element_text(size=8), legend.text=element_text(size = 7))

ggsave("../data/10_clumps_openings_compared.jpg", width = 10, height = 7)

table comparison

current vs target vs prescription

clump_n_trees_grp_summary_temp2 = get_clump_n_trees_grp_summary(
  trees = ttops_temp2, clumps = clump_polys_temp2
)
clump_n_trees_grp_summary_temp2 %>% dplyr::glimpse()
## Rows: 5
## Columns: 31
## $ stand_area_ha               <dbl> 2, 2, 2, 2, 2
## $ clump_n_trees_grp           <ord> Individual, 2-4 trees, 5-9 trees, 10-15 tr…
## $ n_trees                     <int> 64, 39, 56, 21, 57
## $ mean_dbh_cm                 <dbl> 48.15347, 44.26804, 45.02925, 38.01966, 43…
## $ mean_tree_height_m          <dbl> 17.43778, 16.08613, 16.38036, 14.77738, 15…
## $ loreys_height_m             <dbl> 22.07153, 23.74545, 22.43439, 21.57890, 23…
## $ basal_area_m2               <dbl> 14.270633, 8.535523, 11.883499, 3.273830, …
## $ clump_area_ha               <dbl> 0.18087306, 0.08665170, 0.11212985, 0.0341…
## $ stand_n_clumps              <int> 64, 16, 8, 2, 3
## $ trees_per_ha                <dbl> 32.0, 19.5, 28.0, 10.5, 28.5
## $ basal_area_m2_per_ha        <dbl> 7.135316, 4.267761, 5.941749, 1.636915, 5.…
## $ qmd_cm                      <dbl> 53.28276, 52.78832, 51.97963, 44.55264, 51…
## $ stand_trees_per_ha          <dbl> 118.5, 118.5, 118.5, 118.5, 118.5
## $ stand_basal_area_m2         <dbl> 49.95235, 49.95235, 49.95235, 49.95235, 49…
## $ stand_basal_area_m2_per_ha  <dbl> 24.97618, 24.97618, 24.97618, 24.97618, 24…
## $ pct_stand_basal_area        <dbl> 0.28568490, 0.17087328, 0.23789668, 0.0655…
## $ pct_stand_n_trees           <dbl> 0.27004219, 0.16455696, 0.23628692, 0.0886…
## $ stand_qmd_cm                <dbl> 51.80347, 51.80347, 51.80347, 51.80347, 51…
## $ mean_dbh_in                 <dbl> 18.97247, 17.44161, 17.74152, 14.97975, 17…
## $ qmd_in                      <dbl> 20.99341, 20.79860, 20.47998, 17.55374, 20…
## $ stand_qmd_in                <dbl> 20.41057, 20.41057, 20.41057, 20.41057, 20…
## $ mean_tree_height_ft         <dbl> 57.19592, 52.76250, 53.72757, 48.46981, 52…
## $ loreys_height_ft            <dbl> 72.39463, 77.88509, 73.58479, 70.77880, 75…
## $ basal_area_ft2_per_ac       <dbl> 31.102844, 18.603171, 25.900086, 7.135313,…
## $ stand_basal_area_ft2_per_ac <dbl> 108.8712, 108.8712, 108.8712, 108.8712, 10…
## $ trees_per_ac                <dbl> 12.9600, 7.8975, 11.3400, 4.2525, 11.5425
## $ stand_trees_per_ac          <dbl> 47.9925, 47.9925, 47.9925, 47.9925, 47.9925
## $ stand_area_ac               <dbl> 4.942, 4.942, 4.942, 4.942, 4.942
## $ clump_area_ac               <dbl> 0.44693734, 0.21411636, 0.27707286, 0.0842…
## $ basal_area_ft2              <dbl> 153.60909, 91.87636, 127.91398, 35.23951, …
## $ stand_basal_area_ft2        <dbl> 537.6871, 537.6871, 537.6871, 537.6871, 53…
target_data_temp %>% 
  dplyr::select(
    clump_n_trees_grp
    , tidyselect::starts_with("pct_stand_n_trees")
    , tidyselect::starts_with("trees_per_acre_")
  ) %>% 
  dplyr::left_join(
    clump_n_trees_grp_summary_temp2 %>% 
      dplyr::select(clump_n_trees_grp, pct_stand_n_trees, trees_per_ac) %>% 
      dplyr::rename(
        pct_stand_n_trees_presc = pct_stand_n_trees
        , trees_per_acre_presc = trees_per_ac
      )
    , by = "clump_n_trees_grp"
  ) %>% 
  dplyr::mutate(
    dplyr::across(
      tidyselect::where(is.numeric)
      , ~ dplyr::coalesce(.x,0)
    )
  ) %>% 
  dplyr::relocate(
    pct_stand_n_trees_presc
    , .after = pct_stand_n_trees_target  
  ) %>% 
  dplyr::relocate(
    trees_per_acre_presc
    , .after = trees_per_acre_target  
  ) %>% 
  dplyr::mutate(
    dplyr::across(
      tidyselect::starts_with("pct_stand_n_trees")
      , ~ scales::percent(.x,accuracy = 1)
    )
  ) %>%
  kableExtra::kable(
    caption = paste0(
      "Current stand BA (ft2/ac): "
      , clump_n_trees_grp_summary_temp$stand_basal_area_ft2_per_ac[1] %>% scales::number(accuracy = 0.1)
      , "<br>Current stand QMD (in): "
      , clump_n_trees_grp_summary_temp$stand_qmd_in[1] %>% scales::number(accuracy = 0.1)
      , "<br>Current stand TPA: "
      , clump_n_trees_grp_summary_temp$stand_trees_per_ac[1] %>% scales::number(accuracy = 1)
    )
    , escape = F
    , digits = 1
    , col.names = c(
      "", rep(c("current","target","prescription"),2)
    )
  ) %>% 
  kableExtra::kable_styling() %>% 
  kableExtra::add_header_above(
    c(" "=1,"% Trees"=3, "TPA"=3)
  )
Current stand BA (ft2/ac): 221.9
Current stand QMD (in): 16.2
Current stand TPA: 156
% Trees
TPA
current target prescription current target prescription
Individual 2% 30% 27% 2.4 14.4 13.0
2-4 trees 4% 25% 16% 6.1 12.0 7.9
5-9 trees 4% 20% 24% 6.3 9.6 11.3
10-15 trees 1% 10% 9% 2.2 4.8 4.3
16-25 trees 5% 15% 24% 8.1 7.2 11.5
>25 trees 84% 0% 0% 130.6 0.0 0.0

check the pre- and post-treatment summary

# call it
clump_n_trees_grp_summary_temp2 = get_clump_n_trees_grp_summary(
  trees = ttops_temp2
  , clumps = clump_polys_temp2
  )

# table it
clump_n_trees_grp_summary_temp2 %>% 
  dplyr::mutate(presc = "post-treatment") %>% 
  dplyr::bind_rows(
    clump_n_trees_grp_summary_temp %>% 
      dplyr::mutate(presc = "pre-treatment")
  ) %>% 
  dplyr::arrange(clump_n_trees_grp, desc(presc)) %>% 
  dplyr::select(
    clump_n_trees_grp, presc
    , n_trees
    , mean_dbh_in
    , qmd_in
    , mean_tree_height_ft
    , loreys_height_ft
    , trees_per_ac
    , basal_area_ft2_per_ac, pct_stand_basal_area, pct_stand_n_trees
  ) %>% 
  dplyr::mutate(
    dplyr::across(
      .cols = c(pct_stand_basal_area, pct_stand_n_trees) 
      , .fns = ~ scales::percent(.x, accuracy = 1)
    )
  ) %>% 
  kableExtra::kbl(
    digits = 1
    , escape = F
    , caption = "Pre- and Post-treatment overstory tree clump summary"
    , col.names = c(
      ".", ""
      , "trees"
      , "mean<br>DBH (in)"
      , "QMD (in)"
      , "mean<br>Ht. (ft)"
      , "Loreys<br>Ht. (ft)"
      , "TPA"
      , "BA<br>ft<sup>2</sup> ac<sup>-1</sup>"
      , "%<br>stand BA"
      , "%<br>stand trees"
      )
  ) %>% 
  kableExtra::kable_styling() %>% 
  kableExtra::collapse_rows(columns = 1, valign = "top")
Pre- and Post-treatment overstory tree clump summary
. trees mean
DBH (in)
QMD (in) mean
Ht. (ft)
Loreys
Ht. (ft)
TPA BA
ft2 ac-1
%
stand BA
%
stand trees
Individual pre-treatment 12 24.3 26.6 66.2 80.2 2.4 9.3 4% 2%
post-treatment 64 19.0 21.0 57.2 72.4 13.0 31.1 29% 27%
2-4 trees pre-treatment 30 17.8 21.8 52.8 81.4 6.1 15.8 7% 4%
post-treatment 39 17.4 20.8 52.8 77.9 7.9 18.6 17% 16%
5-9 trees pre-treatment 31 13.5 15.7 45.6 65.7 6.3 8.4 4% 4%
post-treatment 56 17.7 20.5 53.7 73.6 11.3 25.9 24% 24%
10-15 trees pre-treatment 11 11.5 12.4 41.9 50.7 2.2 1.9 1% 1%
post-treatment 21 15.0 17.6 48.5 70.8 4.3 7.1 7% 9%
16-25 trees pre-treatment 40 18.3 21.4 54.6 76.7 8.1 20.2 9% 5%
post-treatment 57 17.2 20.4 52.4 75.9 11.5 26.1 24% 24%
>25 trees pre-treatment 645 13.6 15.3 46.2 61.6 130.6 166.3 75% 84%